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Abstract 


The  performance  of  the  dynamic  subgrid-scale  eddy-viscosity  model  and  the  suitabil¬ 
ity  of  high-order  accurate,  upwind-biased  numerical  methods  for  large-eddy  simula¬ 
tions  of  complex  flows  are  investigated  in  the  case  of  the  turbulent  wake  behind  a 
circular  cylinder  at  Reynolds  number  3, 900,  based  on  freestream  velocity  and  cylinder 
diameter. 

The  numerical  method  consists  of  high-order  upwind-biased  finite  difference  tech¬ 
niques  applied  to  the  compressible  Navier-Stokes  equations  written  in  generalized 
coordinates.  Integration  in  time  is  done  using  a  fully  implicit,  second-order  accurate 
iterative  technique.  The  results  of  three  fifth-order  accurate  simulations  performed 
on  identical  grids  with  the  least-squares  version  of  the  dynamic  model,  the  fixed- 
coefficient  Smagorinsky  model,  and  with  no  subgrid-scale  model  are  compared  in  the 
first  10  diameters  of  the  wake.  The  impact  of  three-dimensionality  is  also  examined 
via  two  and  three-dimensional  calculations  without  a  subgrid-scale  model.  The  effect 
of  numerical  dissipation  is  investigated  by  comparing  two  simulations  using  upwind- 
biased  schemes,  the  first  being  fifth-order,  and  the  second  seventh-order  accurate. 

It  is  found  that  the  near-wake  is  highly  three-dimensional  at  this  Reynolds  number. 
It  contains  pairs  of  counter-rotating  streamwise  vortices,  the  effect  of  which  cannot 
be  reproduced  in  two-dimensional  calculations.  Three-dimensional  computations  are 
essential  for  predicting  flow  statistics  of  engineering  interest. 

Amongst  the  three-dimensional  simulations,  although  the  overall  results  are  com¬ 
parable,  the  one  which  uses  the  dynamic  subgrid-scale  model  predicts  more  accurate 
mean  velocities  and  Reynolds  stresses  in  the  vortex-formation  region,  which  includes 
the  first  four  diameters  of  the  wake.  The  fixed-coefficient  model  yields,  overall,  the 


m 


least  accurate  results  in  that  region.  The  magnitude  of  the  eddy-viscosity  in  the  tran¬ 
sitioning  free  shear  layers  and  the  near-wake  is  in  better  agreement  with  expectations 
based  on  the  flow  physics  when  computed  using  the  dynamic  model  rather  than  the 
fixed-coefficient  model.  Differences  between  the  three  calculations  are  largest  in  the 
vortex-formation  zone.  In  the  near-wake  between  4  and  10  diameters  downstream  of 
the  cylinder,  the  extent  and  magnitude  of  the  differences  diminish  due  to  numerical 
dissipation.  Ten  diameters  downstream,  the  three  simulations  predict  comparable 
mean  velocities  and  Reynolds  stresses. 

In  the  near-wake  where  the  mesh  coarsens,  numerical  dissipation  is  found  to  signif¬ 
icantly  affect  the  small  scales  computed  using  the  fifth-order  accurate,  upwind-biased 
scheme.  The  seventh-order  accurate  scheme  predicts  improved  low-order  statistics  at 
the  cylinder  surface  and  in  the  near-wake.  With  the  seventh-order  accurate  scheme, 
smaller-scale  structures  emerge  in  the  wake,  which  are  absent  in  the  simulations  with 
the  fifth-order  scheme.  Velocity  power  spectra  indicate  a  significant  increase  in  energy 
content  at  higher  frequencies  relative  to  the  fifth-order  accurate  calculations,  reflect¬ 
ing  the  impact  of  reduced  levels  of  numerical  dissipation.  Based  on  these  calculations, 
we  have  concluded  that  even  high  order  upwind-biased  schemes  are  overly  dissipative 
and  ill-suited  for  large-eddy  simulations. 
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Chapter  1 
Introduction 


1.1  Review  of  Circular  Cylinder  Flow  Regimes 

For  the  past  century  the  flow  over  a  circular  cylinder  has  been  the  subject  of  numer¬ 
ous  investigations,  of  theoretical,  experimental,  and  more  recently  numerical  variety. 
Extensive  reviews  of  the  knowledge  accumulated  on  this  flow  appear  every  decade 
(Morkovin  1964,  Berger  &  Wille  1972,  Norberg  1987),  yet  our  understanding  of  the 
subjacent  physics  is  incomplete.  The  relevant  non-dimensional  parameter  in  the  flow 
is  the  Reynolds  number,  but  because  of  sensitivity  to  experimental  conditions,  global 
statistics  such  as  drag,  pressure  coefficient,  and  Strouhal  frequency  vary  by  as  much 
as  25%  from  one  experiment  to  the  next  at  a  fixed  Reynolds  number  (Cantwell  & 
Coles  1983).  These  variations  indicate  that  in  addition  to  the  Reynolds  number,  sev¬ 
eral  parameters  arising  from  experimental  set-ups  are  important  in  the  cylinder  flow. 
These  include  the  blockage  ratio,  the  free-stream  turbulence  intensity,  the  cylinder 
aspect  ratio,  and  the  end  boundary  conditions,  each  of  which  has  been  the  subject  of 
numerous  studies. 

A  broad  classification  of  the  cylinder  flow  behavior  in  different  regimes  of  Reynolds 
number  is  presented  in  table  (1),  which  also  lists  some  representative  experiments  for 
each  flow  regime.  The  table  follows  the  descriptions  of  this  flow  presented  by  Morkovin 
(1964),  Roshko  (1954a, b,  1961,  1969),  and  Norberg  (1987). 

At  Reynolds  numbers  less  than  approximately  40,  the  flow  is  laminar  and  steady. 
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The  boundary  layer  on  the  cylinder  surface  separates  at  a  Reynolds  number  of  3.2 
(Nisi  k  Porter  1923)  to  5  (Taneda  1956),  and  a  pair  of  steady  symmetric  vortices  form 
behind  the  cylinder.  A  large  body  of  early  experimental  work  has  documented  that 
range  of  Reynolds  numbers,  detailing  the  main  features  of  the  boundary  layer  and 
near- wake  region  ( e.g .  Tritton  1959,  Coutanceau  k  Bouard  1977,  Thom  1933,  Taneda 
1956,  Acrivos  et  al.  1965).  Between  Reynolds  numbers  10  and  40,  the  velocity  profiles 
in  the  wake  are  self-similar  past  10  diameters  downstream  of  the  cylinder,  the  length 
of  the  recirculation  zone  behind  the  cylinder  grows  linearly  with  Reynolds  number, 
and  the  velocity  distributions  on  the  rear  symmetry  axis  in  the  recirculation  zone  at 
different  Reynolds  numbers  exhibit  similarity  (Nishioka  k  Sato  1974). 

For  Reynolds  numbers  higher  than  approximately  40,  a  characteristic  frequency 
expressed  in  non-dimensional  form  as  the  Strouhal  number  is  associated  with  the 
wake.  Taneda  (1956)  puts  the  critical  Reynolds  number  at  which  shedding  first  occurs 
at  30,  whereas  Homman  (1936),  Kovasznay  (1949)  and  Roshko  (1954b)  find  that 
shedding  starts  at  Reynolds  number  40.  Local  linear  parallel  stability  theory  applied 
to  the  cylinder  wake  (Monkewitz  1988)  indicates  that  the  flow  becomes  absolutely 
unstable  at  Reynolds  number  25,  approximately  two-thirds  the  value  at  which  a 
Karman  vortex  street  develops  experimentally,  showing  that  the  notion  of  absolute 
instability  cannot  by  itself  predict  a  precise  Strouhal  frequency  even  at  the  onset 
of  vortex  shedding.  Chomaz,  Huerre  k  Redekopp  (1988)  use  the  concept  of  global 
instability,  and  show  that  global  oscillations  of  a  shear  flow  will  occur  only  once  a 
critical  sub- volume  of  the  flow  is  absolutely  unstable. 

For  Reynolds  numbers  up  to  150,  the  flow  remains  laminar  (Bloor  1964,  Roshko 
1954a),  the  shed  vorticity  decays  as  it  convects  downstream,  and  the  Strouhal  number 
increases  with  Reynolds  number.  A  least-squares  curve-fit  of  the  Strouhal  curve  for 
Reynolds  numbers  between  50  and  180  is  given  by 

St  =  A/ Re  +  B  +  C  Re  (1) 

where  A  =  -3.3265,  B  =  0.1816  and  C  =  1.6  x  10“4  (Williamson  1991). 

Transition  to  three-dimensionality  in  the  near-wake  occurs  around  Re  —  180,  and 
is  signaled  by  two  discontinuities  in  the  Strouhal- Reynolds  number  relation.  The  first, 
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around  Re  =  180,  arises  from  the  generation  of  vortex  loops  evolving  into  pairs  of 
counter-rotating  streamwise  vortices  in  the  wake,  the  second  comes  from  a  transition 
to  fine-scale  streamwise  vorticity  at  Re  =  230  ~  260. 

At  Reynolds  numbers  between  300  and  2  x  105,  the  sub-critical  range,  the  flow 
around  the  entire  periphery  of  the  cylinder  is  laminar,  and  transition  to  turbulence 
occurs  in  the  separated  free  shear  layers  (Cardell  1993).  At  the  lowest  Reynolds  num¬ 
bers  in  this  range,  the  wake  becomes  fully  turbulent  in  40  to  50  cylinder  diameters 
downstream  (Uberoi  1969),  after  which  distance  the  regular  vortices  have  completely 
decayed.  At  the  higher  end  of  the  Reynolds  number  range,  transition  occurs  very 
near  the  wall  surface,  and  the  wake  is  fully  turbulent  close  downstream  of  the  cylin¬ 
der  (Cantwell  k  Coles  1983).  For  Reynolds  numbers  larger  than  104,  transition  in 
the  shear  layers  occurs  very  close  to  the  separation  points,  and  the  base-pressure  co¬ 
efficient,  drag  coefficient  and  Strouhal  number  are  approximately  constant  at  values 
of  -1.1,  1.2  and  0.2  respectively  (Roshko  k  Fiszdon  1969). 

In  the  past  fifteen  years,  studies  of  the  turbulence  structure  in  plane  wakes  have 
attempted  to  delineate  the  contributions  of  the  organized  and  random  motions  using 
a  triple  decomposition  of  the  flow  into  a  global  mean,  a  phase-averaged  mean  and  a 
random  component  (Perry  k  Lim  1978,  Perry  k  Watmuff  1981,  Hussain  1983,  1986, 
Hayakawa  k  Hussain  1989,  Matsumura  k  Antonia  1993).  Analyses  of  flow  statistics 
using  this  technique  link  the  production  of  random  Reynolds  shear  stress  to  the 
stretching  of  streamwise  vortices  (or  braids)  between  adjacent  span  wise  rollers  (Kiya 
k  Matsumura  1988,  Antonia  et  al.  1987).  Two  representative  experiments  at  low  and 
high  sub-critical  Reynolds  numbers  investigated  the  properties  of  the  coherent  and 
random  Reynolds  stresses  in  the  wake.  Cantwell  k  Coles  (1983)  analyzed  transport 
processes  in  the  near- wake  at  Reynolds  number  1.4  x  105.  Matsumura  k  Antonia 
(1993)  examined  the  intermediate  wake,  from  10  to  40  diameters  downstream,  at 
Reynolds  number  5, 830. 

At  a  Reynolds  number  of  1.4  x  10s ,  Cantwell  k  Coles  found  that  in  the  first  8 
diameters  downstream  of  the  cylinder,  the  coherent  and  random  Reynolds  stresses  are 
of  comparable  magnitude.  The  peak  values  of  the  mean  random  shear  and  streamwise 
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Reynolds  stresses  are  respectively  36%  and  77%  higher  than  their  periodic  counter¬ 
parts.  The  coherent  motion  makes  a  more  important  contribution  to  the  vertical 
velocity  fluctuations,  and  the  periodic  vertical  Reynolds  stress  peaks  are  17%  higher 
than  the  random  component.  As  vortices  are  shed  from  both  sides  of  the  cylinder 
and  turbulent  convective  mixing  occurs  in  the  base  region,  approximately  55%  of  the 
originally  shed  vorticity  is  lost,  and  the  emerging  vortices  convect  downstream  with 
their  centroids  close  to  the  wake  centerline.  In  a  frame  of  reference  moving  with  a 
shed  vortex,  the  flow  can  be  described  as  a  series  of  centers  and  saddle  points.  Near 
these  saddles  lie  peaks  in  turbulent  random  shear  stresses,  and  turbulence  produc¬ 
tion  results  from  the  stretching  of  vorticity  aligned  with  the  separatrices  defining  the 
topology. 

At  the  lower  sub-critical  Reynolds  number  of  5, 830,  Matsumura  &  Antonia  de¬ 
scribe  the  relative  importance  of  the  random  and  periodic  components  of  the  flow.  As 
in  the  case  of  Cantwell  &  Coles  at  higher  Reynolds  number,  the  coherent  motion  con¬ 
tributes  a  large  portion  of  the  vertical  Reynolds  stresses.  Ten  diameters  downstream, 
random  stresses  account  for  only  25%  of  the  total  vertical  Reynolds  stress,  but  for 
70%  of  the  streamwise  Reynolds  stress.  The  Reynolds  shear  stress  at  the  same  lo¬ 
cation  is  closely  equi-partitioned  between  the  periodic  and  random  components,  and 
its  peak  value  occurs  in  the  centers,  not  in  the  saddle  points  as  is  found  at  Reynolds 
number  1.4  x  105. 

Much  of  the  literature  dealing  with  experiments  on  the  cylinder  documents  phe¬ 
nomena  in  the  sub-critical  Reynolds  number  region.  It  has  been  known  since  the  work 
of  Gerrard  in  1965  that  the  mean  aerodynamic  properties  of  the  cylinder  are  particu¬ 
larly  sensitive  to  free-stream  disturbances  in  the  sub-critical  range  of  Reynolds  num¬ 
ber.  These  disturbances,  as  well  as  acoustic  noise  levels,  cylinder  vibrations,  surface 
roughness,  blockage  ratio,  and  other  geometric  parameters  which  have  not  been  care¬ 
fully  documented  in  most  experiments,  have  been  shown  by  Norberg  (1987),  amongst 
others,  to  influence  transition  to  turbulence  in  the  free  shear  layer,  as  well  as  mixing 
and  entrainment  in  the  wake  region  for  sub-critical  and  critical  Reynolds  numbers  up 
to  3  x  105.  This  may  account  for  the  lack  of  agreement  in  the  literature  on  the  near¬ 
wake  properties  of  the  flow.  The  influence  of  such  perturbations  is  best  illustrated  by 
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the  dependence  of  mean  cylinder-surface  statistics  on  surface  roughness  for  Reynolds 
numbers  larger  than  10s  (Shih,  Wang,  Coles  k  Roshko  1992).  Surface  roughness  is 
quantified  by  a  parameter  representing  the  mean  protrusion  height  ( k )  at  the  cylinder 
surface,  normalized  by  the  cylinder  diameter  ( D ).  At  a  Reynolds  number  of  106,  the 
mean  drag  and  base-pressure  coefficients  of  smooth  cylinders  are  approximately  0.25 
and  —0.3  respectively.  For  k/D  =  0.0003,  the  mean  drag  increases  by  a  factor  of  3 
to  0.75,  and  the  base-pressure  coefficient  drops  to  —0.8.  As  (k/D)  further  increases 
to  0.01,  the  drag  rises  to  1.1,  while  the  base-pressure  coefficient  falls  to  —1.3.  Ex¬ 
perimental  results  of  Achenbach  (1971),  Roshko  (1961),  and  Shih,  Wang,  Coles  k 
Roshko  (1992)  suggest  that  in  the  extreme  case  of  very  rough  cylinders,  some  mean 
statistics,  including  the  base-pressure  coefficient,  the  pressure  rise  near  separation, 
and  the  boundary  layer  separation  angle,  become  Reynolds-number  independent  and 
vary  only  with  surface  roughness  beyond  a  Reynolds  number  of  2  x  106. 

The  critical  range  of  Reynolds  numbers,  between  2  x  10s  and  3.5  x  106,  displays  two 
transitions  in  the  drag  coefficient,  labeled  the  lower  and  upper  transitions  by  Roshko 
(1961).  In  the  lower  transition  range  (2  x  105  <  Re  <  5  x  105),  the  drag  coefficient 
drops  abruptly  from  1.2  to  about  0.3  due  to  an  increase  in  base  pressure  at  a  Reynolds 
number  of  approximately  3.6  x  105.  A  laminar  separation  of  the  boundary  layer  is 
followed  by  transition  to  turbulence,  reattachment  and  a  final  turbulent  separation. 
The  separation  point  moves  from  the  front  to  the  downstream  side  of  the  cylinder, 
and  the  width  of  the  near-wake  decreases  to  less  than  1  diameter.  In  the  upper 
transition  region  (5  x  10s  <  Re  <  3.5  x  106),  the  base-pressure  coefficient  decreases 
monotonically  from  approximately  —0.2  to  —0.5,  while  the  drag  coefficient  increases 
from  0.3  to  0.7,  and  remains  at  that  value  for  Reynolds  numbers  of  up  to  107.  With 
increasing  Reynolds  number  the  separation  point  moves  forward,  but  it  remains  on 
the  downstream  side  of  the  cylinder,  and  the  wake  width  increases  but  stays  smaller 
than  1  diameter.  The  sensitivity  of  the  flow  to  disturbances  in  the  critical  regime, 
in  particular  the  non-zero  mean  lift  which  can  develop  around  the  cylinder,  has  been 
experimentally  investigated  by  Schewe  (1986). 

In  the  post-critical  regime,  past  Reynolds  number  3.5  x  106,  the  boundary  layer 
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on  the  cylinder  surface  transitions  to  turbulence  before  separating.  The  separation- 
reattachment  bubble  present  in  the  critical  region  disappears.  The  base-pressure 
coefficient  pursues  its  monotonic  decrease  started  at  Reynolds  number  5  x  105,  reach¬ 
ing  -0.6  at  Reynolds  number  8  x  106;  the  drag  coefficient  is  constant  at  around  0.7, 
and  vortices  are  shed  regularly  at  an  approximately  constant  Strouhal  frequency  of 
0.27  (Roshko  1961). 

1.2  On  Comparison  with  Experiments  and  the 
Case  Selected  for  Simulations 

The  flow  in  the  sub-critical  Reynolds  number  range  (300  <  Re  <  2  x  105)  features 
several  challenging  phenomena  for  large-eddy  simulations:  a  laminar  boundary  layer 
including  unsteady  separations  and  reattachments,  flow  reversals  at  the  cylinder  sur¬ 
face  and  in  the  near-wake,  adverse  pressure  gradients,  transitioning  free  shear  layers, 
and  a  turbulent  wake  with  random  and  periodic  Reynolds  stresses  of  comparable 
magnitudes. 

The  selection  of  a  particular  Reynolds  number  represents  a  compromise  between 
competing  effects.  At  higher  Reynolds  numbers,  random  motions  account  for  an  in¬ 
creasingly  larger  share  of  the  total  Reynolds  stresses  close  to  the  cylinder,  but  the 
laminar  cylinder  boundary  layer  and  the  transitioning  separated  shear  layers  become 
thinner.  The  resolution  requirements  of  these  layers  increase  with  the  Reynolds  num¬ 
ber,  independently  of  whether  the  turbulent  wake  is  simulated  via  direct  or  large-eddy 
simulations. 

The  choice  of  Reynolds  number  is  also  dependent  on  the  available  experimental 
data.  Accurate  measurements  of  velocities  and  stresses  are  difficult  in  low  or  reversed 
mean  velocity  regions  where  turbulence  intensities  are  high,  as  well  as  in  high  flow- 
angle  regions,  which  accounts  for  the  dearth  of  complete,  published  experimental  data 
in  the  near  wake.  A  careful  review  of  the  literature  reveals  that  despite  the  imposing 
volume  of  published  work  on  the  circular  cylinder,  most  experiments  concentrate 
on  the  mean  and  fluctuating  properties  at  the  cylinder  surface  ( e.g  Achenbach  1968, 
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Roshko  1954a,  West  1990,  Farell  k  Blessmann  1983,  Schewe  1983,  Norberg  1987),  the 
intermediate  wake  region  between  10  and  40  diameters  downstream  (e.g.  Matsumura 
k  Antonia  1993,  Zhou  k  Antonia  1992,  Hayakawa  k  Hussain  1989),  or  the  small- 
deficit  far-wake  region  extending  from  50  to  several  hundred  diameters  downstream 
( e.g  Townsend  1949,  Antonia  k  Browne  1987,  Wygnanski  et  a 1.  1986,  Fabris  1979, 
Papailiou  k  Lykoudis  1974,  Ferre  et  a 1.  1990). 

Experiments  detailing  near-wake  properties  at  sub-critical  Reynolds  numbers  are 
few,  and  all  but  one  provide  inadequate  or  insufficient  data  for  comparison  with  a 
numerical  simulation.  Bouard  k  Coutanceau  (1980)  study  the  time  evolution  of  the 
near- wake  velocity  profiles  at  Reynolds  numbers  between  40  and  104,  but  provide  no 
Reynolds-stress  information.  Bloor  (1964),  Wei  k  Smith  (1986)  and  Kourta  et  a 1. 
(1987)  concentrate  on  transitional  waves  in  the  separated  shear  layers  and  do  not 
discuss  other  statistics.  Bloor  k  Gerrard  (1966)  provide  mean  streamwise  velocity 
profiles,  at  Reynolds  numbers  of  2, 000  and  16, 000,  at  different  locations  downstream 
of  the  recirculation  region,  but  give  no  details  on  the  Reynolds  stresses,  the  mean 
vertical  velocities,  or  the  velocity  distribution  inside  the  bubble. 

One  experiment  which  provides  detailed  Reynolds  stress  and  velocity  information 
within  the  stagnation  zone  and  in  the  near- wake  region  was  performed  by  Cantwell 
and  Coles  (1983),  at  a  Reynolds  number  of  1.4  x  105.  Mean  and  phase-averaged 
velocity  and  Reynolds  stress  fields  are  documented  within  8  diameters  downstream 
of  the  cylinder,  as  well  as  mean  and  fluctuating  quantities  at  the  cylinder  surface. 
At  the  lower  Reynolds  number  of  3,900,  two  separate  experiments  provide  details  on 
the  velocity  and  Reynolds-stress  distributions.  Lourenco  et  al.  (1993)  used  a  Particle 
Image  Velocimetry  technique  to  compile  mean  and  phase- averaged  data  within  three 
diameters  downstream  of  the  cylinder.  In  the  near-wake  region,  between  the  closure 
point  of  the  recirculation  bubble  and  10  diameters  downstream,  single  sensor  mea¬ 
surements  of  the  mean  velocities  and  Reynolds  stresses  have  been  documented  by  Ong 
and  Wallace  (1994).  The  data  at  both  Reynolds  numbers  of  3,900  and  1.4  x  105  are 
adequate  for  comparison  with  large-eddy  simulation  results.  On  a  structured  mesh, 
and  with  the  fifth-order  spatially-accurate  numerical  method  used  in  this  work,  the 
number  of  mesh  points  in  the  plane  of  mean  motion  required  at  the  higher  Reynolds 
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number  is  approximately  one  order  of  magnitude  larger  than  at  Reynolds  number 
3, 900,  a  point  further  discussed  in  chapter  4.  For  this  reason,  the  lower  Reynolds 
number  is  selected  for  comparison  with  large-eddy  simulations. 

1.3  On  Computing  the  Flow  Over  Circular  Cylin¬ 
ders 

The  many  physical  phenomena  arising  in  the  turbulent  flow  around  circular  cylinders 
(section  1.1)  represent  a  difficult  challenge  for  numerical  prediction  methods.  At  the 
1980-1981  HTTM-Stanford  Conference  on  Complex  Turbulent  Flows,  the  experiment 
of  Cantwell  k  Coles  (1983),  documenting  mean  flow  quantities  in  the  near  wake  at 
Reynolds  number  1.4  x  105,  was  suggested  as  a  test  case  for  numerical  methods.  Of 
the  52  participating  computational  groups,  none  attempted  simulating  this  case,  for 
reasons  which  were  not  documented.  In  recent  years,  the  advent  of  more  powerful 
computers  and  sophisticated  numerical  methods  has  led  to  new  efforts  to  simulate 
the  transitional  and  turbulent  cylinder  flow  regimes.  The  features  and  principal 
conclusions  of  several  representative  calculations  are  discussed  below. 

In  the  transitional  Reynolds  number  range  (150  <  Re  <  300),  two-dimensional 
simulations  can  accurately  predict  a  restricted  number  of  global  parameters,  such  as 
mean  drag  on  the  cylinder  and  Strouhal  number,  without  modeling  three-dimensional 
effects.  Braza,  Chassaing  k  Ha  Minh  (1986)  use  a  finite-volume  method,  second- 
order  accurate  in  space  and  time,  to  discretize  the  two-dimensional  Navier-Stokes 
equations  written  in  logarithmic-polar  coordinates.  At  Reynolds  numbers  of  100  and 
200,  the  computed  Strouhal  frequency  and  mean  drag  coefficient  are  in  good  agree¬ 
ment  with  experimental  measurements.  Franke,  Rodi  k  Schonung  (1990)  attempt 
to  predict  these  coefficients  at  Reynolds  numbers  as  high  as  5,000  with  a  spatially 
third-order  accurate,  first-order  time-accurate,  finite-volume  method  applied  to  the 
equations  written  in  standard  polar  coordinates.  Having  confirmed  the  grid  and 
time-step  independence  of  their  results,  they  noted  that  Strouhal  frequency  and  drag 
coefficient  could  be  accurately  computed  only  for  Reynolds  numbers  under  300.  Past 
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this  Reynolds  number,  in  the  sub-critical  regime,  the  computed  drag  increases  with 
Reynolds  number,  contrary  to  the  experimental  trend,  which  the  authors  surmise  is 
a  consequence  of  neglecting  the  three-dimensionality  of  the  flow. 

Visualization  experiments  (Bays-Muchmore  &  Ahmed  1992,  Williamson  1989, 
1991b)  demonstrate  the  presence  of  counter-rotating  streamwise  vortices  in  the  braid 
region  between  consecutive  spanwise  vortices,  directly  behind  the  cylinder,  for  Reynolds 
numbers  higher  than  approximately  200.  The  importance  of  accounting  for  the  con¬ 
sequent  three-dimensionality  is  further  corroborated  by  the  results  of  Karniadakis  & 
Triantafyllou  (1992),  who  performed  direct  numerical  simulations  of  the  cylinder  flow 
at  Reynolds  numbers  between  175  and  500  with  a  mixed  spectral/spectral-element 
technique.  The  time-trace  of  the  velocity  at  a  representative  point  in  the  near-wake, 
calculated  from  a  two-dimensional  simulation  at  Reynolds  number  500,  exhibits  a 
quasi-periodic  behavior.  In  the  three-dimensional  simulations  however,  the  flow  be¬ 
havior  alternates  between  apparently  random  and  almost  periodic  states,  which  is  in 
qualitative  agreement  with  the  experimental  observations  of  Bloor  (1964). 

These  low  Reynolds  number  simulations  investigated  the  onset  of  three-dimensional 
motions  in  the  wake  as  well  as  the  mechanisms  involved  in  transition  to  turbulence. 
Three-dimensionality  in  the  wake  was  found  to  arise  from  a  secondary  instability  of 
the  vortex  street  appearing  in  the  computations  between  Reynolds  numbers  of  200 
and  210,  and  was  not  accompanied  by  the  hysteresis  phenomenon  documented  ex¬ 
perimentally  by  Williamson  (1988).  Transition  to  turbulence  is  analyzed  from  power 
spectra  and  three-dimensional  phase  portraits  at  locations  in  the  near-wake.  These 
indicate  that  a  period-doubling  bifurcation  cascade  may  be  responsible  for  transition 
to  turbulence.  However,  time-averaged  statistics  at  the  cylinder  surface  and  in  the 
near- wake  are  not  compared  with  experimental  data,  making  it  difficult  to  assess  the 
accuracy  of  the  computations. 

Steady  Reynolds-averaged  computations  at  higher  than  transitional  Reynolds 
numbers  were  undertaken  by  Majumdar  &  Rodi  (1985),  who  used  a  k  —  e  model 
to  compute  sub-critical  and  post-critical  cases,  at  Reynolds  numbers  of  1.4  x  105  and 
3.6  x  106  respectively.  Majumdar  &  Rodi  sought  to  establish  the  level  of  inaccuracy 
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of  predictions  based  on  steady  RANS  simulations.  Their  method  consisted  of  a  finite- 
volume  formulation,  with  a  hybrid  central/upwind  scheme  applied  to  the  convective 
terms.  Neither  the  mean  forces  acting  on  the  cylinder,  nor  the  near-wake  mean  flow 
were  predicted  accurately.  At  the  sub-critical  Reynolds  number  of  1.4  x  105,  the 
solution  was  compared  with  the  experiments  of  Achenbach  (1968)  and  Cantwell  k 
Coles  (1983).  The  quantities  which  compared  well  with  these  experiments  are  the 
wall  shear  stress  upstream  of  separation,  and  the  wake  spreading-rate  downstream  of 
the  stagnation  region.  However  most  other  relevant  quantities  were  predicted  incor¬ 
rectly.  The  boundary-layer  separation  occured  too  far  downstream,  the  base  pressure 
was  too  high  and  the  drag  coefficient  correspondingly  too  low.  In  the  wake  region, 
the  maximum  shear  stress  was  4  times  smaller,  and  the  recirculation  bubble  80% 
longer  than  found  experimentally.  These  results  were  qualitatively  insensitive  to  the 
treatment  of  the  eddy  viscosity  in  the  laminar  and  attached  part  of  the  cylinder 
boundary  layer.  Setting  the  eddy  viscosity  to  zero  in  that  region,  or  using  a  Van 
Driest  damping  function  resulted  in  an  improved  wall  shear  stress  distribution  past 
the  separation  point.  At  the  post-critical  Reynolds  number  of  3.6  x  106,  separation 
was  predicted  11%  too  far  downstream  along  the  cylinder  surface,  the  base  pressure 
was  too  high  although  the  minimum  pressure  coefficient  compared  fairly  well  with 
experimental  measurements,  and  the  wall  shear  stress  was  far  too  large,  with  a  peak 
2.8  times  higher  than  found  experimentally.  Similar  discrepancies  were  reported  by 
Sugawanam  k  Wu  (1982). 

Seeking  the  cause  of  the  difficulties  encountered  in  these  computations,  Franke, 
Rodi  k  Schonung  (1989)  studied  the  phase- averaged  experimental  data  of  Cantwell  k 
Coles  at  Reynolds  number  1.4  x  105,  in  which  they  examined  the  validity  of  the  k  —  e 
formulation.  Their  principal  finding  was  that  the  eddy  viscosity  is  anisotropic,  and 
is  negative  in  regions  where  history  and  transport  effects  dominate  over  production 
of  Reynolds  stresses,  indicating  that  in  the  case  of  the  circular  cylinder,  the  k  —  e 
model  is  not  adequate.  Algebraic  stress  models  have  been  applied  to  flows  around 
square  cylinders  (Franke,  Rodi  k  Schonung  1989,  Franke  k  Rodi  1991,  Murakami  et 
a 1.  1991,  Murakami  1990),  but  not  to  circular  cylinders. 
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The  difficulties  encountered  by  Reynolds-averaged  models  have  recently  moti¬ 
vated  the  undertaking  of  three-dimensional  simulations  of  the  sub-critical  and  critical 
regimes.  These  efforts  have  been  restricted  to  demonstrating  improvements  in  the 
predictions  of  the  mean  and  fluctuating  forces  acting  on  the  cylinder  when  three- 
dimensionality  is  introduced  in  the  simulations,  and  do  not  concentrate  on  the  sta¬ 
tistical  description  of  the  wake  region.  At  Reynolds  number  104,  a  finite-element 
large-eddy  simulation  was  performed  by  Kato  &  Ikegawa  (1991),  using  the  standard 
Smagorinsky  model  and  a  linear  wall-damping  function.  The  computed  mean  drag 
coefficient  (1.14),  base-pressure  coefficient  (-1.08),  Strouhal  number  (0.20)  and  fluc¬ 
tuating  lift  coefficient  (0.27)  are  in  good  agreement  with  experimental  values.  In  a 
similar  effort  to  demonstrate  the  improvements  in  these  coefficients  when  extracted 
from  three-dimensional  simulations,  Tamura,  Ohta  &  Kuwahara  (1990)  performed 
finite-difference  calculations  at  Reynolds  numbers  of  104,  105  and  106  without  tur¬ 
bulence  models.  The  incompressible  Navier-Stokes  equations,  written  in  generalized 
coordinates  in  the  plane  of  mean  motion,  are  discretized  with  an  upwind,  third-order 
accurate  scheme  for  the  convective  terms,  and  marched  in  time  with  a  semi-implicit 
first-order  accurate  method.  The  mean  drag  is  in  reasonable  agreement  with  ex¬ 
perimental  results  at  the  aforementioned  Reynolds  numbers,  although  the  computed 
base-pressure  coefficient,  of  —0.5  at  Reynolds  number  106,  is  far  lower  than  the  ex¬ 
perimental  value  of  —0.3. 

The  collective  findings  of  the  simulations  described  in  this  section  are  summa¬ 
rized  in  table  (2).  Two-dimensional,  unsteady  Navier-Stokes  simulations  can  predict 
Strouhal  numbers  and  drag  coefficients  at  transitional  Reynold  numbers,  but  be¬ 
come  unreliable  in  the  sub-critical  regime.  Direct  numerical  simulations,  performed 
to  examine  the  physics  of  the  transition  to  three-dimensionality  and  turbulence  in 
the  near-wake,  have  been  limited  to  maximum  Reynolds  numbers  of  500.  Steady 
Reynolds-averaged  (RANS)  simulations  using  the  k  —  e  model  at  sub-critical  and 
post-critical  Reynolds  numbers  predict  inaccurate  velocity  and  Reynolds  stress  dis¬ 
tributions  in  the  near-wake,  while  giving  mixed  results  at  the  cylinder  surface.  Coarse 
direct  simulations  improve  on  the  RANS  calculations,  predicting  reasonable  drag  co¬ 
efficients  up  to  a  critical  Reynolds  number  of  106,  providing  further  evidence  that 


11 


resolving  the  three-dimensionality  of  the  flow  is  of  foremost  importance  past  the 
transitional  regime.  A  large-eddy  simulation  at  the  low  sub-critical  Reynolds  number 
of  104  was  able  to  predict  both  drag  and  wall- pressure  coefficients  in  agreement  with 
experiments,  although  an  accurate  resolution  of  the  near-wake  was  not  attempted. 


1.4  Present  Status  of  Large-Eddy  Simulation 

In  large-eddy  simulation  the  dynamically  important  large-scale  motions  are  computed 
directly,  while  a  subgrid-scale  turbulence  model  represents  the  effect  of  the  unresolved 
small  scales  on  the  large  scales. 

In  the  Smagorinsky  subgrid-scale  eddy-viscosity  model,  the  eddy  viscosity  is  lo¬ 
cally  proportional  to  the  large-scale  strain-rate  tensor,  and  includes  an  adjustable 
coefficient.  A  brief  review  of  the  adjustments  required  in  different  flows  is  presented 
by  Germano,  Piomelli,  Moin  k  Cabot  (1991),  who  note  that  although  specific  compli¬ 
cating  factors  in  a  given  flow  may  be  accommodated  by  altering  the  model  constant, 
no  single  coefficient  performs  well  under  arbitrary  flow  conditions.  Further  empiricism 
is  introduced  into  the  Smagorinsky  model  by  the  use  of  wall  damping  and  intermit- 
tency  functions,  which  ensure  the  proper  asymptotic  behavior  of  the  subgrid-scale 
stresses  near  solid  boundaries  and  in  transitioning  flows. 

The  mismatch  between  the  premise  that  small  scales  tend  to  be  more  universal 
than  large  ones,  and  the  degree  of  arbitrariness  prevalent  in  early  LES  applications 
motivated  the  use  of  modern  statistical  turbulence  theories  in  the  search  for  more 
general  models.  A  description  of  the  models  which  subsequently  emerged  in  the  1980’s 
and  early  1990’s  is  given  by  Moin  k  Jimenez  (1993).  The  model  of  Yakhot  k  Orszag 
(1986)  based  on  renormalization  group  theory  has  incorrect  limiting  wall  behavior;  in 
a  modified  version  developed  by  Piomelli  et  a 1.  (1990),  grid-independent  results  are 
found  difficult  to  achieve.  Chollet  k  Lesieur  (1981),  following  Kraichnan’s  (1976)  eddy 
viscosity  model  in  spectral  space,  proposed  a  simplified  model  formulated  in  Fourier 
space.  This  model  cannot  be  used  in  complex  flow  applications.  A  related  model, 
developed  by  Metais  k  Lesieur  (1992)  to  remedy  this  shortcoming  and  expressed  in 
physical  space,  does  not  have  the  correct  asymptotic  behavior  near  walls  and  is  nearly 
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identical  to  Smagorinsky’s  model. 

The  concept  of  dynamic  subgrid-scale  modeling,  introduced  by  Germano,  Piomelli, 
Moin  &  Cabot  in  1991,  is  based  on  the  scale-similarity  ideas  of  Bardina  et  al.  (1980), 
in  which  spectral  information  available  from  the  resolved  scales  is  used  to  model 
the  small  scales.  When  used  in  conjunction  with  the  Smagorinsky  formulation,  the 
Smagorinsky  coefficient  becomes  a  function  of  space  and  time,  dynamically  computed 
from  the  resolved  large-scale  fields.  The  distinguishing  properties  of  the  dynamic 
eddy  viscosity  model  are  its  correct  asymptotic  behaviors  near  solid  surfaces,  and  its 
intrinsic  ability  to  differentiate  between  laminar  and  turbulent  regions  of  the  flow, 
without  a d  hoc  damping  or  intermittency  functions. 

The  dynamic  subgrid-scale  model  has  been  generalized  to  compressible  flows  and 
scalar  transport  by  Moin,  Squires,  Cabot  &  Lee  (1991)  and  Cabot  &;  Moin  (1991).  The 
method  for  computing  the  coefficient  of  the  model,  initially  restricted  in  applicability 
to  flows  statistically  homogeneous  in  at  least  one  direction,  has  been  generalized  by 
Ghosal,  Lund  &  Moin  (1992)  to  apply  to  inhomogeneous  flows  using  a  constrained 
optimization  formulation.  This  model,  tested  in  the  decay  of  homogeneous  isotropic 
turbulence  and  in  channel  flow,  gives  results  in  good  agreement  with  the  experimental 
data. 

Various  formulations  of  the  dynamic  model  have  been  successfully  tested  against 
experimental  and  direct  numerical  simulation  results  in  a  variety  of  flows  involv¬ 
ing  specific  complications.  These  include  transitional  and  fully  developed  turbulent 
channel  flow  (Germano,  Piomelli,  Moin  &  Cabot  1991),  channel  flow  with  passive 
scalars  (Cabot  &  Moin  1991),  isotropic  turbulence  decay  (Moin,  Squires,  Cabot  & 
Lee  1991,  Ghosal,  Lund  &  Moin  1992),  planetary  boundary  layer  (Bohnert  &  Ferziger 
1993),  three-dimensional  channel  flow  (Cabot  1993),  channel  flow  with  system  rota¬ 
tion  (Squires  &  Piomelli  1993  and  Cabot  1993),  boundary  layers  with  imbedded 
streamwise  vortices  (Liu  &  Piomelli  1993),  flow  over  a  backward-facing  step  (Ak- 
selvoll  &  Moin  1993a, b)  and  driven  cavity  flows  (Zang,  Street  h  Koseff  1993).  The 
quality  of  the  results  achieved  in  these  simulations  suggests  that  the  dynamic  pro¬ 
cedure  is  highly  successful  for  the  simulation  of  flows  in  which  a  variety  of  physical 
phenomena  are  present  in  different  sub-regions. 
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1.5  Objectives  and  Overview 


The  first  objective  of  this  study  is  to  evaluate  the  performance  of  the  dynamic  subgrid- 
scale  eddy  model  in  a  complex  flow  where  Reynolds-averaged  turbulence  models  have 
faced  difficulties.  One  of  the  most  challenging  applications  of  dynamic  models  to  date 
has  been  the  computation  of  the  flow  over  a  backward-facing  step  by  Akselvoll  k  Moin 
(1993a,b).  The  turbulent  wake  behind  a  circular  cylinder  is  the  first  application  of 
dynamic  models  to  external  flows.  It  involves  unsteady  separation  and  reattachment, 
laminar  sub-regions  and  transition  to  turbulence.  Past  attempts  at  simulating  the 
sub-critical  regime  with  a  k  -  t  model  reveal  the  inability  of  the  Reynolds-averaged 
approach  to  accurately  predict  the  near- wake  statistics.  The  previous  results  of  large- 
eddy  simulations  using  the  Smagorinsky  model  and  coarse  direct  simulations  were 
restricted  in  scope  to  evaluating  the  forces  acting  on  the  cylinder.  The  sub-critical 
cylinder  wake  flow  is  then  appropriately  suited  for  a  validation  study  of  a  dynamic 
subgrid-scale  eddy  viscosity  model,  which  adapts  to  local  flow  conditions  without 
using  adjustable  constants. 

Assessing  the  importance  of  three-dimensional  effects  is  another  objective  of  this 
work.  To  this  end,  the  near-wake  and  cylinder-surface  statistics  obtained  from  two- 
dimensional  and  large-eddy  simulations  are  compared. 

The  contribution  of  the  subgrid-scale  eddy  viscosity  model  per  se  is  evaluated  by 
comparing  mean  velocities  and  Reynolds  stresses  computed  from  simulations  without 
a  turbulence  model,  with  the  dynamic  subgrid-scale  model,  and  with  a  fixed- coefficient 
Smagorinsky  model. 

Finally,  one  of  our  objectives  is  to  evaluate  the  suitability  of  high-order  accu¬ 
rate,  upwind-biased  finite-difference  methods  for  the  large-eddy  simulation  of  flows  in 
complex  geometry.  Results  of  two  simulations  which  use  the  dynamic  subgrid-scale 
turbulence  model  are  compared  to  evaluate  the  impact  of  upwinding  on  the  computed 
solution.  One  uses  a  fifth-order  accurate,  one  point  up  wind- biased  method  (Rai  k 
Moin  1991),  the  other  a  seventh-order  accurate,  one  point  upwind-biased  scheme  to 
discretize  the  convective  derivatives.  In  both  cases,  viscous  fluxes  are  evaluated  using 
sixth-order  accurate,  central  differencing  schemes. 
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This  study  describes  one  of  several  concurrent  efforts  at  Stanford  to  evaluate 
various  formulations  of  dynamic  models  in  complex  flows.  These  include  a  curved 
boundary  layer  at  Reynolds  number  1,200  based  on  momentum  thickness,  the  flow 
around  an  airfoil  at  Reynolds  number  106,  and  the  flow  through  a  plane  diffuser  at 
Reynolds  number  9, 000.  For  further  information  on  these  calculations,  the  reader  is 
referred  to  forthcoming  reports  published  by  the  Center  for  Turbulence  Research. 

The  findings  and  contributions  of  this  report  are  as  follows: 

•  A  high-order  accurate,  upwind-biased  numerical  method  was  developed  for  solv¬ 
ing  compressible  flow  problems  in  generalized  coordinates. 

•  At  Reynolds  number  3,900,  the  near- wake  behind  the  circular  cylinder  con¬ 
tains  pairs  of  counter-rotating  streamwise  vortices  between  the  spanwise  rollers.  It  is 
strongly  three-dimensional,  and  cannot  be  accurately  simulated  with  a  two-dimensional 
calculation. 

•  At  the  same  Reynolds  number,  a  large-eddy  simulation  using  the  least-squares 
version  of  the  dynamic  subgrid-scale  eddy-viscosity  model  was  overall  more  accurate 
in  the  first  four  diameters  of  the  wake  than  its  counterparts  using  the  fixed-coefficient 
Smagorinsky  model  or  no  subgrid- scale  model. 

•  The  magnitude  of  the  eddy  viscosity  in  the  separated  free  shear  layers  and  the 
near-wake  region  is  in  better  agreement  with  expectations  based  on  the  flow  physics 
when  using  the  dynamic  model  rather  than  the  fixed-coefficient  model. 

•  The  numerical  dissipation  generated  by  the  fifth-order  accurate,  upwind-biased 
scheme  applied  to  the  convective  fluxes,  was  found  to  have  a  significant  impact  on 
the  turbulence  in  the  near- wake. 

•  The  simulation  using  the  seventh-order  accurate,  upwind-biased  scheme  applied 
to  the  convective  terms  demonstrated  the  effects  of  decreased  numerical  dissipation 
on  the  computed  solution.  Smaller-scale  structures  appeared  in  the  near-wake  that 
were  not  present  when  using  the  fifth-order  accurate  upwind-biased  scheme,  and  the 
energy  of  the  computed  solution  was  substantially  increased  at  higher  frequencies. 


1.6  Summary 


The  next  chapter  describes  the  numerical  method  used  in  all  the  simulations  in  this 
report.  Chapter  3  presents  the  test-problems  validating  the  method  and  its  imple¬ 
mentation:  the  laminar,  two-dimensional  flow  over  a  circular  cylinder  at  steady  and 
unsteady  Reynolds  numbers,  as  well  as  the  linear-stability  analysis  of  a  forced  channel 
flow  in  three  dimensions.  Large-eddy  simulations  of  the  flow  over  a  circular  cylinder 
at  Reynolds  number  3, 900  are  documented  in  chapter  4,  where  the  dynamic  eddy 
viscosity  model  is  evaluated  and  compared  to  the  Smagorinsky  model.  Comparisons 
of  the  results  obtained  with  the  fifth  and  seventh-order  accurate,  upwind  biased  dif¬ 
ferencing  schemes  are  presented  in  chapter  5.  Conclusions  regarding  the  performance 
of  the  dynamic  model  and  the  suitability  of  the  numerical  method  are  presented  in 
chapter  6,  where  recommendations  for  extensions  of  this  work  are  provided. 
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Table  1:  Summary  of  cylinder  flow  regimes 
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Turbulent  cylinder  boundary  layer  Roshko  1961 

Post-critical  >  3.5  x  106  Regular  vortex  shedding  (St  ~  0.27),  Co  ~  0.7  Shih  et  a/.  1992 

Transition  precedes  separation,  no  reattachment 


Re  Range  — > 


l  Simulation  J, 


Two-dimensional 

Unsteady 

Navier-Stokes 


RANS 

with 

k  —  t  model 


Transitional 

regime 


Cj),  St 


Sub-critical 

regime 


Unreliable 


Wall  shear  ahead 
of  separation 


Critical  I  Post-critical 


regime 


regime 


Wall  Cp  up 
to  recovery 


Inaccurate  near-wake  velocity 
and  Reynolds  stress  distributions 


CD,  St,  Wall  CP 
(Wake  not  resolved, 
Re  =  104) 


Onset  of  3-D 


Transition 
(Re  =  500) 


CD 

(Wake  not 
resolved) 


Table  2:  Status  of  prediction  capabilities  for  flow  over  a  cylinder 
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Chapter  2 

Mathematical  Formulation 


2.1  Large-Scale  Equations  of  Motion 


The  equations  governing  large-scale  motions  of  the  flow  are  derived  by  first  considering 
the  continuity,  momentum  and  total  energy  equations: 


d  d  ,  .  A 

a?’W'>'“)  =  0 

d  d 

-Qt{pui)  +  ^  +  pSik  -  <?ik )  =  o 

d  d 

~die  +  ^xl^Uk^e  +P^~ <TikU'  +  =  0 

where  the  molecular  heat  flux  and  stress  tensor  are  given  by 

7  dT 


Qi  = 


(7  —  1  )RePr  dx{ 


1  (dui  duk\  2  1  duj 
o%k  =  —  I  - - V  ~  I  -  ~—<>ik - 


;) 


Re\dxk  '  dx{ )  3  Re'  dx/ 


(2) 

(3) 

(4) 

(5) 

(6) 


The  fluid  is  assumed  to  obey  the  perfect-gas  law.  The  total  energy  and  pressure  are 
thus  related  through 

7}  1 

e  = - -  +  -pukuk  (7) 

7  —  1  z 


These  equations  are  normalized  with  free-stream  sound-speed  Coo  and  cylinder  radius 
Rc.  The  quantity  Re  refers  to  the  Reynolds  number,  Re  =  p<x>c<x,Rc/ PI  the  molecular 
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viscosity  p  is  assumed  to  be  independent  of  temperature;  the  Prandtl  number  Pr  and 
Mach  number  are  fixed  at  0.72  and  0.2  respectively  in  all  cylinder  simulations. 

The  compressible  formulation  is  chosen  to  facilitate  the  potential  inclusion  of 
zonal  methods  into  the  numerical  scheme.  The  principal  issue  in  these  methods 
is  the  prescription  of  boundary  conditions  at  zonal  interfaces.  Boundary  conditions 
based  on  momentum  and  energy  fluxes  across  sub- volumes  can  be  easily  implemented 
for  compressible  methods,  and  have  been  successfully  applied  to  a  variety  of  flow 
situations  (Rai  1985,  1986a, b). 

The  large-scale  equations  are  obtained  by  applying  a  filtering  operation,  denoted 
by  an  overbar,  to  the  governing  equations.  A  filtered  quantity  <f>  is  defined  as: 


m  =  [  G(x,mv)dy 

JD 


where  D  is  the  entire  flow  domain.  G  is  a  spatial  filter  which  removes  high  spatial- 
frequency  information,  with  characteristic  length  in  direction  x,  given  by  A,-.  Filtering 
the  equations  leads  to: 

^+d^m  =  0  (9) 

d  d 

— pu7  -(-  - — (pUkUi  +  pf>ik  ~  Vik)  =  0.  (10) 

Ot  OXk 

y  +  h  (”‘<e  +  ^  -  (7  -  Sf)  =  (11) 

Momentum  fluxes  are  decomposed  into  resolvable  and  modeled  components  as 

_  WiWk  ,  T _  puiWk]  no\ 

PUkUi  =  — n - b  pUkUi - z —  (1^) 

pi  pi 


where  Tjfc  is  the  subgrid-scale  stress  tensor.  Following  Moin  et  a 1.  (1991),  the  subgrid- 
scale  viscous  work  is  assumed  to  be  negligible  in  the  energy  balance,  leading  to: 


GikV'i  —  &ik  V* i 


The  subgrid-scale  transport  term  to  be  modeled  in  the  filtered  energy  equation: 

Uk{e  +  P)  -ujb(e-fp), 
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upon  using  (7)  for  the  total  energy,  can  be  written  as: 

_  >*y  ^ 

uk(e  +  p)  ~  «fc(e  +  p)  = - r qk  +  -(puj^-Uj  -  UkpUjUj) 

—  1  Z 


where  qk  is  the  subgrid-scale  heat  flux: 

qk=W^-^  =  W~k-lJ^L  (16) 

P 

The  second  term  on  the  right-hand  side  of  (15)  represents  the  subgrid-scale  transport 
of  total  kinetic  energy,  which  is  smaller  than  the  subgrid-scale  heat  flux  by  a  factor 
proportional  to  M It  is  neglected  for  the  present  low  Mach  number  applications. 
The  large-scale  equations  of  motion  are  then: 

|?+^=°  (17) 


d _  ,  d  ( pu{  puk 

dt^U%  dxk  \  p 


+  pSik  -  (Tik  +  Tik  )  =  0 


JL-.JL 

dt€  dxk 


(ut(e + T) -  w  -  (7  _  + 7^T®)  =  0  (19) 


2.2  Turbulence  Models 

The  large-scale  momentum  and  energy  equations  are  closed  by  modeling  the  subgrid- 
scale  Reynolds  stress  tensor  r,j  and  heat  flux  vector  qk-  Comparisons  will  be  presented 
in  Chapter  4  between  closures  obtained  from  dynamic  modeling  and  the  Smagorinsky 
model  with  a  damping  function  and  a  fixed  coefficient.  Each  approach  is  described 
in  the  remainder  of  this  Section. 

In  the  first  approach,  subgrid-scale  terms  are  evaluated  with  the  local  least-squares 
version  of  dynamic  turbulence  models  (Moin  et  al.  1991,  Lilly  1992),  in  which  model 
coefficients  are  only  functions  of  time  and  azimuthal  and  radial  directions.  These 
models  use  large-scale  information  extracted  from  the  computed  fields  via  a  test¬ 
filtering  operation.  In  the  following  derivations,  a  test-filtered  function  /  is  denoted 
by  the  symbol  /.  In  this  work,  the  test-filter  is  chosen  as  a  tophat  filter  in  physical 
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space,  and  is  applied  in  all  spatial  directions.  The  trapezoidal  rule  is  used  to  evaluate 
the  resulting  volume  integrals  numerically.  The  characteristic  scales  of  the  grid  and 
test-filtering  operations  are  denoted  by  A  and  A  respectively. 


2.2.1  Subgrid-Scale  Reynolds  Stresses 

A  subgrid-scale  stress  tensor  Tij  at  the  test-filter  scale  A  is  defined  by  analogy  with 
the  subgrid-scale  Reynolds  stress  r,j  as 


m  Pui  Pui 

_Z —  pu^uj  ^ 

P 


(20) 


The  subgrid-scale  stress  tensors  at  both  grid  and  test-filter  levels  are  modeled  using 
a  trace-free  Smagorinsky  formulation: 

ry-^r»  =  -2C(*,!/,t);sAJ|3|35  (21) 

Tij  -  f  =  -2 C(x,  y, tff  A2jfjf|  (22) 

where  S*j  is  the  trace-free  rate  of  strain  tensor, 

=  (23) 


and  |5|  is  the  norm  of  Sif  _ 

Pi  =  y/tSifri.  (24) 

C(x,  t)  is  a  model  coefficient  to  be  determined.  T,j  and  f,j  can  be  related  through  the 
Germano  identity  to  the  Leonard  stress  tensor  Lif 


pui  pUj 

7 


(25) 


where  L,j  is  computable  from  the  resolved  large-scale  field.  The  model  coefficient 
C(x,  y,  t)  is  determined  by  substituting  (21)  and  (22)  into  (25).  The  resulting  relation 

Lij  =  2C(x,y,t)A2  Mij  (26) 
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where 


(27) 


is  solved  for  C  by  contracting  both  sides  with  the  tensor  M,j,  which  is  equivalent  to 
a  least-squares  minimization  procedure  (Lilly  1992).  Denoting  the  spatial  average  in 
the  spanwise  direction  by  <  >  ,  CA2  is  given  by: 


C(x,y,t)  A2 


1  <  LjjMjj  > 

2  <  MijMji  >' 


(28) 


For  reasons  of  numerical  stability,  negative  values  of  the  model  coefficient  C  are 
discarded  by  an  ad  hoc  clipping  operation,  which  results  in  a  modified  coefficient  C *: 

C*(x,y,t)A2  =  -(|CA2|  +  CA2)  (29) 


The  only  adjustable  constant  in  the  model  is  the  ratio  of  the  test  to  grid-filter 
widths  A/A.  In  all  calculations,  this  ratio  is  chosen  to  be  2  (Germano  et  a 1991). 
The  turbulent  eddy  viscosity  is  found  from  (21)  as: 


fit  =  C*A2p  |S|. 


(30) 


2.2.2  Subgrid-Scale  Energy 

In  compressible  flow,  the  isotropic  part  of  the  subgrid-scale  Reynolds  stress  tensor 
rkk  =  q2  cannot  be  combined  with  the  thermodynamic  pressure.  Following  Moin  et 
a 1.  (1991),  q2  is  modeled  using  the  relation: 

q 2  =  2Ci(x,  i i,  t)p  A2|*Sj2.  (31) 

Contracting  (25)  and  substituting  (31)  and  a  similar  expression  for  Tkk  leads  to  an 
expression  for 


Averaging  this  expression  in  the  spanwise  direction  and  substituting  in  (31)  gives  the 
subgrid-scale  energy  q 2: 


The  expressions  (29)  and  (33)  define  the  subgrid-scale  Reynolds  stress  tensor  ry, 
closing  the  filtered  momentum  equations. 


2.2.3  Subgrid-Scale  Heat  Flux 

A  model  for  subgrid-scale  heat  flux  qk  is  needed  to  close  the  energy  equation.  An 
eddy-diffusivity  model  is  used,  which  introduces  a  turbulent  Prandtl  number  Prt  to 
be  determined  dynamically: 


4fc  =  ~P 


vt  dT  _  _C*(x,  y,  t)A2|6'|  dT 
Prt  dxk  ^  Prt  dxk 


Writing  the  heat  flux  at  the  test-filter  level  as  Hk  using  the  same  eddy-diffusivity 


model: 


rr  ~  fM  -C*A2|5|  dT 
Hk-puk  %  ~  P  Pr  dxh' 


the  difference  Hk  —  qk  is  computable  from  the  large-scale  field: 


Hk-qk  = 


ppuk  p  puk 


C*A2  ((A\*  _^dT  ^OT 

tytU a) 


Following  Moin  et  a  1.  (1991),  the  expression  for  the  turbulent  Prandtl  number  is  found 
by  contracting  the  equation  above  with  dT/dxk  and  averaging  in  the  homogeneous 
spanwise  direction: 


Prt  =  C*(x,y,t)  A 


<  |P0jUfc 


p  J  &xk 
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This  expression  closes  the  large-scale  energy  equation.  The  filtered  equations 
and  associated  subgrid-scale  models  have  been  presented  in  physical  space.  The 
corresponding  discrete  equations  to  be  solved  numerically  are  formulated  in  compu¬ 
tational  space.  The  transformation  between  the  physical  and  computational  spaces 
is  presented  in  Section  2.3. 

2.2.4  Fixed-Coefficient  Smagorinsky  Model 

In  simulations  which  use  a  fixed-coefficient  Smagorinsky  model,  the  subgrid-scale 
energy  q 2  (equation  31)  is  neglected  (i.e.  t**  =  0)  and  the  turbulent  Prandtl  number 
is  set  to  1.  The  residual  stress  is: 

Tij  =  -2 C,p  A2 1  S*  I  SjT.  (38) 

The  eddy  viscosity  fit  is  related  to  the  length  scale  /  of  the  model  by: 

l>,  =  C.A>  |  S7  |=  Pp  |  S=  |  .  (39) 

The  length  scale,  chosen  to  match  that  used  by  Piomelli,  Ferziger  &  Moin  (1987),  is 
given  in  generalized  coordinates  by: 

l  =  Ca\fl  -  e-(r~1)+3/>‘+3  (J  A(AqAz)1^3  (40) 

with  Cs  =  0.065  and  A+  =  25.  The  Jacobian  J  of  the  two-dimensional  coordinate 
transformation  from  ( x,y )  to  (^,?/)-space  is  defined  in  (42),  r  =  r((,q)  denotes  the 
radial  distance  from  the  wall,  £  is  the  wall  normal  direction.  Wall  coordinates  are 
defined  using  the  mean  skin-friction  coefficient  C j  to  provide  a  velocity  scale: 

(41) 

The  Mach  number  appears  in  this  relation  since  the  Reynolds  number  is  defined  with 
the  freestream  sound  speed  as  a  velocity  scale. 


K  -  1)+  =  «  -  l)MooRe 
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2.3  Transformation  to  Generalized  Coordinates 

The  governing  equations  in  physical  space  (x,y,z,t)  are  transformed  to  their  ho- 
mologues  in  computational  space  (£,  17,  z,  t)  via  a  two-dimensional  time-independent 
mapping.  In  the  cylinder  simulations,  £,  9  and  z  correspond  to  generalized  radial, 
azimuthal  and  span  wise  directions  respectively.  Letting  J  represent  the  Jacobian  of 
this  coordinate  transformation 

J  —  •Pjj  3/(5  (42) 


derivatives  in  physical  and  computational  spaces  are  related  by  the  expression 


(d/dx\  /  y„  -yA  ( d/d£\ 

\d/dy)  '  \-xn  x£  )  \a/dv) 


(43) 


The  computational  space  is  a  cube  with  equispaced  point  distributions  in  all  direc¬ 
tions.  The  mesh  spacing  in  the  transformed  £  and  9  directions  are  set  to  =  A9  =  1 
for  convenience,  while  the  spacing  in  the  spanwise  direction  is  determined  by  the  span- 
wise  box  and  mesh  sizes.  Denoting  as  Q  the  conservative-variable  vector 


f  P  > 

92 

pu 

93 

=  J 

~pv 

94 

pw 

V95/ 

\  e  / 

(44) 


the  large-scale  equations  can  be  formally  written  as 


d_ 

at 


~  a  ^ 

o+§ff«+ 


^  p  ,  ^  p  ~  0 

&jF- + &f*  - 0 


where  F £  represents  the  total  flux  along  the  (  direction.  Detailed  expressions  for  the 
fluxes  appearing  in  (45)  are  given  in  Appendix  A. 
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2.4  Spatial  Discretization 

2.4.1  Flux-Splitting  and  Non-Conservative  Viscous-Flux  For¬ 
mulation 

The  governing  equations  in  generalized  coordinates  (45)  were  presented  in  terms  of 
the  total  fluxes  F^,  Fn  and  Fz  along  each  direction.  Each  of  these  fluxes  can  be  written 
as  the  sum  of  convective,  viscous  and  heat  flux  components.  In  a  generic  direction  x, 
these  are  given  by  F®,  F”  and  Fx  respectively: 

Fx  =  F®  +  F”  +  Fx  (46) 

Before  differencing  these  three  fluxes  numerically,  it  is  necessary  to  rewrite  each  into 
different  forms  as  follows. 

The  inviscid  ’Euler’  flux  F®  is  decomposed  into  signals  propagating  in  opposite 
directions  along  all  axes  of  the  computational  grid  using  flux-vector  splitting  (Steger 
&  Warming  1981).  The  flux  is  diagonalized  as 

Fx®  =  PX\XP~1Q.  (47) 

where  Ax  represents  the  diagonal  eigenvalue  matrix.  The  total  convective  flux  is  then 
written  as 

Fex  =  F*+  +  FT  =  Px\iP;'Q  +  PxKP;lQ  (48) 

where 

A*  =  j(A,  ±  |A,|).  (49) 

Viscous  and  heat  fluxes  are  decomposed  into  three  components  each,  Fjx,  Fnx 
and  Fzx,  which  contain  terms  involving  only  primitive  variable  derivatives  along  the 
t)  and  z  directions  respectively.  Thus  in  each  flux  vector  FX,X2,  the  first  index 
(xj)  indicates  that  the  vector  contains  derivatives  of  primitive  variables  along  the  xx 
direction  only,  while  the  second  index  (#2)  indicates  the  direction  of  the  flux.  For  the 
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two-dimensional  coordinate  mapping  described  in  the  previous  section,  the  viscous 
and  heat  fluxes  are: 


II 

jS1 

+ 

J'qi 

+ 

& 

II 

✓ty  ^ 

+ 

-**31 

(50) 

f;  =  4 + 4 + 4;  4  =  4 + 4 

(51) 
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For  numerical  stability  reasons,  the  second  derivative  operators  acting  on  the  primitive 
variables  must  be  expressed  separately  from  first  derivative  operators.  Second  deriva¬ 
tives  appear  in  the  governing  equations  through  the  viscous  and  heat  fluxes,  which 
are  functions  of  Q  and  its  spatial  derivatives:  F#  =  F^(Q,Q^),  Fm  —  Fv,t(Q,  Qv) 
and  Fzz  =  FZZ(Q,  Qz)-  Each  of  these  fluxes  can  be  written  as 

F,,(QA)  =  aJ4  (53) 


where  Axx  =  dFxx/dQx.  The  flux  Fxx  and  is  operated  on  in  (45)  by  the  operator 
d/dx,  leading  to 

dFxx  .  d2Q  dAxxdQ  , , .  x 

~fa~Axxdx*  +  dx  dx  •  1  ’ 


This  chain-rule  expansion  yields  non-conservative  viscous  flux  terms. 

Combining  these  transformations  in  the  expressions  for  the  convective,  viscous 
and  heat  fluxes,  the  governing  equations  (45)  become: 


upwind  inviscid  fluxes 


downwind  inviscid  fluxes 


In  ,  j?e+  ,  JL  pe+  ,  d_pe+  ,  jL pc-  ,  jLpe~  +  A/?e- 

dtQ  +  di Fi  +dv  '  +  dz  z  +dt  *  +dv  '  +  dz  z 

+|(4  +  +  pi)  +  +  4)^ 

+^(4 +K + 4) + (K, + <0 + ^ 

+  Tzfa + + (A-* + +  i(A" + A‘*)  &  = 
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Each  flux  and  matrix  appearing  in  (55)  is  described  in  Appendix  A.  Numerical 
discretization  of  this  equation  requires  the  use  of  two  distinct  spatial  differencing 
schemes  to  be  applied  to  the  upwind  and  downwind  convective  fluxes,  as  well  as  first 
and  second  derivative  schemes  for  the  various  components  of  the  viscous  terms. 

2.4.2  Convective  and  Viscous  Differencing  Schemes 

In  the  present  work,  convective  derivatives  are  evaluated  with  two  high-order  accu¬ 
rate,  one-point  upwind-biased  schemes  in  two  different  simulations.  The  fifth-order 
accurate  scheme  (Rai  &  Moin  1991)  is  given  on  a  uniform  mesh  at  grid  location  x,  by 

^  f  I  ■‘■r  i  *■  r  t  i  *  t  i  *  t  *  t  i  “  ■>  I  (5g) 

(57) 
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1  r  1  r  r  1  f  ,  1  f  1  ,  ,  A 6  d6/ 

“  ~3(r‘~3  +  4r’~2  ~  •'‘-1  +  3^‘  +  2*i+1  ~  20 •'i+2  +  60  dx& 

1  1 1  l  r  l  r  l  A6  d6f 

i  =  20^~2  ~  2^,_1  ~  3^  +  ^+1  ~~  4^*+2  +  30 ^,+3  _  60"  dx* 


where  A  represents  the  mesh  spacing.  The  seventh-order  accurate  scheme  is 


d+f 


dx+ 
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The  first  and  second  derivatives  appearing  in  the  viscous  terms  are  evaluated  using 
central,  sixth-order  accurate  schemes: 
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=  ™ (/«'-3  +  /*'+3)  -  ™(/»- 2  +  /i+2)  +  +  /»+i)  ~  Tq/* 


A8 


(60) 


(61) 


90  w 20v,,,_‘  '  '  2W*“‘  ‘  •’“ri/  18J‘  560  dxs 

Since  the  schemes  applied  to  the  convective  terms  are  biased  one  point  upwind,  in¬ 
troducing  sixth-order  accurate  first  derivatives  does  not  increase  the  global  stencil 


size. 
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2.4.3  Modified  Wave-Numbers 

Applying  the  convective  upwind  schemes  to  a  function  f(k)e'kx,  letting  0  =  &A,  and 
denoting  the  modified  wave-number  of  the  scheme  by  0  =  &A,  the  fifth  and  seventh- 
order  accurate  schemes  yield  0  respectively  given  by: 

0  =  ^sin30-^sin20  +  ^sin0-i(-^cos30-|-^cos20-icos0-|--)  (62) 

>  =  -  T5o sin  4<>  +  xls sin  _  ! sin 

+  i  ^  cos  £0  +  —  cos  30  —  —  cos  20  +  —  cos  0  —  — )  (63) 

Real  and  imaginary  parts  of  these  modified  wave-numbers  are  displayed  in  figures 
(1)  and  (2)  respectively.  The  imaginary  parts  of  the  wave-numbers  arises  because 
of  the  upwind  biasing.  Their  presence  introduces  a  damping  in  an  otherwise  purely 
oscillatory  system.  This  is  clearly  illustrated  by  considering  the  one-dimensional 
wave  equation  du/dt  +  du/dx  =  0.  A  Fourier  mode  of  u  can  be  written  as  f(t)e'kx, 
where  f(t)  =  f(0)e~ikt.  Substituting  for  k  the  modified  wave-number  of  the  scheme 
%  =  kr  +  iki  introduces  in  the  response  of  f(t)  the  term  ek,t.  For  the  fifth-order 
scheme,  since  k  is  negative  at  wave  numbers  larger  than  0.357r/A  (figure  2),  where 
the  differencing  scheme  starts  being  inaccurate,  this  exponential  term  damps  the 
response  at  frequencies  which  are  too  high  to  be  accurately  differentiated  by  the 
scheme  on  a  given  mesh.  The  real  part  of  0  is  compared  in  figure  (1)  with  the 
modified  wave-numbers  of  the  fourth-order  central  scheme: 

/'  =  -  8/j-.  +  S/i+1  -  /j-i)  (M) 


which  is 


6  =  -(8sin0  —  sin  20), 
6 


as  well  as  with  the  modified  wave-number  of  the  fourth-order  accurate  Pade  scheme 

fj+ 1  +  fi-i  +  =  A^j+1  _  ^-1  ^  +  30  j 
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given  by 


(67) 


~  3sin0 
2  +  cos  0 

The  fourth-order  Pade  and  fifth-order  upwind  schemes  have  nearly  identical  modified 
wave-numbers  up  to  k  ~  0.457r/A.  At  higher  wave-numbers  the  error  in  the  Pade 
scheme  is  consistently  smaller  than  that  in  the  upwind  scheme.  The  Pade  scheme  has 
a  purely  real  modified  wave  number,  indicating  that  the  scheme  does  not  generate 
numerical  dissipation  on  a  uniform  grid. 

The  modified  wave- number  for  the  second  derivative  scheme  (61)  applied  to  vis¬ 
cous  fluxes 

io  aq 

d2  =  —  —  cos  30  +  —  cos  20  —  3  cos  0  +  —  (68) 

45  10  18 

is  shown  in  figure  (3)  compared  with  the  modified  wave-number  of  a  second-order  ac¬ 
curate  central  difference  scheme,  for  which  02  =  2(1  -  cos  0).  The  second-order  scheme 
is  accurate  for  wave-numbers  below  approximately  0.3tt/A,  whereas  the  sixth-order 
scheme  is  close  to  the  exact  solution  for  wave-numbers  up  to  0.6x/A.  Wave-numbers 
higher  than  these  limiting  values  are  under-damped  in  the  respective  schemes,  but 
the  sixth-order  method  provides  50%  more  dissipation  than  the  second-order  one  at 
the  highest  wave-number  of  k  =  tt/A. 

Near  the  wall  and  outer  boundaries,  spatial  differencing  schemes  must  be  modified 
with  a  corresponding  reduction  in  order  of  accuracy.  At  these  boundaries,  convective 
and  viscous  derivatives  are  second-order  accurate.  Details  of  the  schemes  used  in 
these  regions  are  provided  in  Appendix  (A). 


2.5  Time  Integration 


The  governing  equations  are  marched  forward  in  time  using  a  variable  time-step,  fully 
implicit  scheme  defined  at  time-level  n  +  1  by: 


dQn+1  a  +  /? 


Qn+1  + 


dt 


a/3 


P 


■Qn  + 


a 


a(a  -  /3)  P(/3  -  o) 


Qn- 


(69) 


where  a  =  Atn+1  is  the  time  step  at  time  level  n  +  1,  and  /3  =  Atn+1  +  At".  For 
constant  time-steps  this  scheme  is  second-order  accurate.  In  the  turbulent  cylinder 
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calculations  presented  in  Chapter  4  and  Appendix  F,  time  steps  were  defined  by 
constant  Courant  numbers.  The  mean  variation  in  At  from  step  to  step  was  ap¬ 
proximately  0.01%  for  a  typical  simulation.  The  time  derivative  dQn+1  /dt  appearing 
in  (69)  is  given  by  the  governing  equations  (55).  Equation  (69)  represents  a  sys¬ 
tem  of  coupled  non-linear  equations.  It  is  solved  at  each  time  step  with  a  modified 
Newton-Raphson  iteration  technique  (Rai  &  Moin  1993).  An  iterative  approach  was 
chosen  because  the  direct  construction  and  inversion  of  the  implicit  systems  defined 
by  (69)  present  several  difficulties.  The  fifth-order  upwind  spatial  differencing  scheme 
leads  to  block  hepta-diagonal  implicit  matrices,  which  are  expensive  to  generate  and 
invert.  In  three  dimensions,  the  size  of  the  implicit  systems  requires  the  use  of  a 
factored  approach.  However  a  direct  factored  scheme  involves  two  additional  difficul¬ 
ties  which  can  restrict  the  time-step  size.  The  viscous  fluxes  which  involve  spatial 
cross-derivatives  (e.g.  d2Q  /  d^dif)  must  be  evaluated  explicitly  (Beam  &  Warming 
1978),  generating  a  time  error  of  order  dependent  on  the  type  of  explicit  scheme,  and 
the  factoring  error,  of  order  At2,  cannot  be  eliminated. 

An  iterative  technique  can  circumvent  these  difficulties.  The  iterative  scheme  is 
based  on  expressing  the  solution  vector  field  Qn+1  at  time-level  tn+1  as  the  zero  of 
the  function 

f[Qn+1]  =  -Qn+1  +  ClQn  +  c2Qn~'  +  c3jtQn+1  (70) 

which  is  a  statement  equivalent  to  (69).  The  coefficients  c\  through  C3  represent  the 
corresponding  combinations  of  a  and  0  in  (69).  The  root  of  /  is  sought  using  a 
classical  Newton  approach: 

-  (4/[«])V+1  -  Q')=  /[<?'].  p!) 


where  p  is  the  iteration  index.  The  Jacobian  on  the  left-hand  side  is  obtained  from 
(70)  as 

(72) 

The  full  iterative  scheme  is  written  by  substituting  (72)  and  (70)  into  (71),  leading 


to: 


(73) 
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where  the  right-hand  side  vector  is 

Rn’p  =  -Qp  +  cxQn  +  c2Qn~l  +  cZjtQp.  (74) 

The  time  derivative  of  Q  in  this  relation  is  given  by  (55).  Upon  convergence  of  the 
iterations,  Qv+l  =  Qp  =  Qn+1 ,  and  the  left-hand  side  of  (73)  is  zero  regardless  of  any 
approximation  made  in  evaluating  the  Jacobian  terms  d(dQ/dt)/dQ.  The  implicit 
system  (73)  is  inverted  by  factoring  the  equations.  The  factoring  error  is  also  driven 
to  zero  upon  convergence  of  the  iterations.  Schematically,  the  three  resulting  factored 
systems  are  written  as: 

(I  +  c3Mv)pQ*  =  Rn*  (75) 

(/  +  c3Mz)pQ ~  =  Q*  (76) 

(/  +  c3MQpAQp+1  =  Q**  (77) 


where  Mv,  Mz  and  represent  linearized  Jacobians  in  the  77,  z  and  £  directions 
respectively.  The  crux  of  the  method  lies  in  approximating  these  Jacobians  so  that 
they  are  inexpensive  to  compute,  define  easily  invertible  matrices,  and  lead  to  stable 
iterations.  The  construction  of  the  Jacobians  in  this  work  follows  the  discussions  of 
Rai  k  Moin  (1991)  and  Pulliam  k  Chaussee  (1981).  At  each  iteration,  viscous  stresses 
are  treated  implicitly  in  the  wall-normal  direction  £,  but  are  neglected  in  the  periodic 
azimuthal  (77)  and  spanwise  (z)  directions.  The  iterative  scheme  finally  becomes: 

(7 + + <**;%)’ >  =  p;'8”  (re) 

(7 + c3a+£ + c3a;| j)  V - p;'p/r  (79) 
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(80) 


The  eigenvector  matrices  Pv  and  obtained  from  the  diagonalization  of  the  az¬ 
imuthal  and  spanwise  Euler  fluxes,  are  described  in  Appendix  A.  On  the  left-hand 
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side  of  these  systems,  convective  derivatives  are  approximated  with  upwind  first-order 
accurate  schemes,  while  central  second-order  accurate  schemes  are  applied  for  viscous 
derivatives.  This  results  in  tridiagonal  scalar  systems  in  the  periodic  directions  and 
a  block-tridiagonal  system  in  the  wall-normal  direction,  which  are  inverted  directly. 

In  all  laminar  and  turbulent  computations  in  this  report,  the  number  of  iterations 
per  time-step  is  set  to  three,  which  is  shown  to  be  sufficient  to  reduce  the  residual  norm 
(|/?"’p|)  by  a  factor  of  102  to  104  at  each  time  step  (Appendix  D,  Section  D.3.3  and 
Appendix  F,  Section  F.3.2).  Validation  tests  in  the  laminar  vortex-shedding  regime 
of  the  flow  over  a  circular  cylinder  (Chapter  3)  provide  some  indication  that  this  level 
of  convergence  is  sufficient  to  ensure  that  the  residual  error  does  not  noticeably  affect 
the  accuracy  of  the  computed  solution. 

Applied  to  the  model  equation 

^  =  A  u  (81) 

at 


the  time  advancement  scheme  (69)  is  stable  for  a  given  time-step  At  provided  the 
inequality 

|  i2  _  ]2_±vTT2AAf£  <  l  (82) 

11  |3  —  2XAt\2  ~  V  ' 

holds.  The  relation  above  also  defines  the  stability  criterion  for  the  model  equation 


du  du 

di+cai=0 


(83) 


applied  to  one  Fourier  component  f(k)e'kx  of  u,  provided  the  coefficient  A  is  set  to 
the  complex  expression 

A  =  -^(^sin30-^sin20+^sin0-i(--^cos30+-cos20--cos0+-)^  (34) 

where  the  term  between  parentheses  represents  the  modified  wave  number  of  the  fifth- 
order  accurate  upwind  convective  scheme  (62).  The  envelope  (\rAxAt/c,  \tAxAt/ c) 
corresponding  to  |<r|2  =  1  is  displayed  in  figure  (4).  The  scheme  is  unconditionally 
stable  for  Ar  <  0. 
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2.6  Boundary  Conditions 


The  boundary  conditions  needed  for  the  simulation  of  the  flow  over  a  circular  cylinder 
are  wall  and  far-field  boundary  conditions  applicable  to  the  conservative  variable 
increment  A Q  (see  equation  80).  Spanwise  and  azimuthal  directions  being  periodic, 
the  appropriate  boundary  conditions  on  the  intermediate  fields  Q*  (equation  78)  and 
Q**  (equation  79)  are  periodic  conditions.  Wall  and  far-field  conditions  are  imposed 
implicitly  for  improved  numerical  stability,  and  are  updated  after  each  time-step  to 
prevent  drifting.  Their  precise  formulation  is  presented  in  Appendix  B. 

At  the  cylinder  surface,  no-slip  boundary  conditions  are  applied.  The  pressure 
derivative  is  evaluated  from  the  normal-momentum  equation,  and  the  surface  heat- 
flux  is  set  to  zero. 

Free  surfaces  are  divided  into  potential  and  wake-region  boundaries,  each  one  of 
which  admits  a  different  set  of  boundary  conditions.  The  conditions  imposed  in  the 
potential  region  are  based  on  locally  one-dimensional  Riemann  invariants  (Pulliam 
1986b)  and  neglect  all  viscous  contributions.  At  the  subsonic  outflow,  the  value  of 
the  first  Riemann  invariant  is  imposed  as  that  given  by  the  potential  flow  around  a 
circular  cylinder.  At  the  subsonic  inflow,  the  required  four  boundary  conditions  on 
entropy,  spanwise  and  tangential  velocities,  and  a  Riemann  invariant,  are  similarly 
based  on  the  potential  flow  solution. 

The  edges  of  the  wake  layer  at  the  domain  boundary  are  determined  using  a 
699  boundary-layer  criterion.  Within  the  wake  outflow  region,  primitive  variables  are 
extrapolated  with  a  first-order  accurate  scheme.  These  boundary  conditions  are  found 
to  be  stable,  but  are  inaccurate  and  generate  perturbation  waves  which  propagate 
from  the  outer  boundary  toward  the  center  of  the  computational  domain.  In  all 
simulations,  a  long  region  of  highly  stretched  mesh  minimizes  the  influence  of  the 
boundary  conditions  on  the  region  of  interest. 


2.7  Numerical  Dissipation  and  Aliasing  Control 


The  upwind-biased  scheme  chosen  to  evaluate  the  derivative  of  the  convective  fluxes 
in  the  equations  of  motions  introduces  numerical  dissipation  in  the  computation.  This 
dissipation  can  affect  the  accuracy  of  the  computed  solution  in  regions  where  it  is  of 
the  same  order  as  the  sum  of  its  molecular  and  subgrid-scale  counterparts.  On  the 
other  hand  numerical  dissipation  is  instrumental  in  stabilizing  the  numerical  method, 
by  dissipating  the  energy  content  of  the  high  frequencies  in  the  flow,  and  thereby 
acting  as  a  built-in  aliasing  control  (not  removal)  mechanism. 

In  direct  simulations  of  incompressible  turbulent  channel  flow  at  Reynolds  number 
180  based  on  wall  shear-velocity  and  channel  half-width,  Rai  k  Moin  (1991)  observed 
that  the  numerical  dissipation  introduced  by  the  upwind-biased  scheme  presented 
above  (equation  56),  was  apparently  large  enough  to  control  aliasing  by  dissipating 
the  high  frequency  content  of  the  solution  fields.  However  a  precise  analysis  of  numer¬ 
ical  dissipation  was  not  provided.  Further  the  ability  of  a  scheme  to  control  aliasing 
errors  in  an  incompressible  calculations  may  not  extend  to  compressible  calculations, 
the  latter  involving  many  more  divisions  and  multiplications  of  the  primitive  variables 
than  the  former.  It  was  observed  from  the  vorticity  fluctuations  that  small-scale  mo¬ 
tions  were  damped  in  a  coarse  simulation  which  resolved  the  mean  flow  accurately. 
Energy  spectra  for  this  coarse  simulation  were  not  presented.  Root  mean  square 
vorticity  fluctuations  agreed  well  with  spectral  calculation  results  in  a  fine-mesh  sim¬ 
ulation.  The  one-dimensional  streamwise  energy  spectra  in  that  simulation  showed 
that  the  energy  at  high  frequencies  was  about  two  and  a  half  decades  below  the  val¬ 
ues  obtained  from  spectral  calculations.  The  spectral  simulation  used  for  comparison 
required  (192  x  129  x  160)  grid-points,  while  the  upwind-biased  finite-difference  cal¬ 
culations  used  a  (192  x  101  x  192)-point  mesh.  Thus  although  the  study  by  Rai  k 
Moin  concluded  that  high-order  upwind-biased  finite-difference  methods  were  good 
candidates  for  direct  simulations  of  complex-geometry  flows,  the  applicability  of  such 
methods  to  large-eddy  simulations  was  not  investigated. 

In  a  large-eddy  simulation  the  prime  objective  is  to  accurately  resolve  low-order 
statistics  on  a  coarse  mesh.  The  use  of  an  upwind-biased  scheme  constrains  the  grid  to 
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be  fine  enough  that  the  large  scales  are  not  affected  by  numerical  dissipation.  Further 
in  such  a  simulation,  a  turbulence  model  will  impact  the  computed  large  scales  only 
if  the  magnitude  of  the  numerical  dissipation  is  sufficiently  small  compared  to  that 
of  the  molecular  dissipation. 

In  the  simulations  at  Reynolds  number  3,900  the  numerical  dissipation  gener¬ 
ated  by  the  fifth-order  upwind-biased  scheme  was  found  to  significantly  damp  the 
one- dimensional  frequency  spectra  between  five  and  ten  diameters  downstream  of  the 
cylinder.  The  level  of  numerical  dissipation  present  in  the  computed  solution  over¬ 
whelmed  the  turbulence  model,  which  did  not  have  an  impact  on  the  mean  velocities 
and  Reynolds  shear  stresses  in  that  region.  The  simulation  with  the  seventh-order 
accurate  upwind-biased  scheme  demonstrated  that  reducing  numerical  dissipation  led 
to  more  energetic  small  scales  in  the  near-wake.  These  results  axe  further  discussed 
in  chapters  4  and  5. 


2.8  Triple  Flow  Decomposition 

The  flow  in  the  cylinder  near-wake  is  characterized  by  a  periodic  vortex  shedding 
motion,  which  at  sub-critical  Reynolds  numbers  contributes  between  a  half  and  a 
third  of  the  total  Reynolds  stresses  (Cantwell  &  Coles  1983,  Matsumura  &  Antonia 
1993).  To  differentiate  between  the  random  and  periodic  components  of  the  Reynolds 
stresses,  a  flow  variable  can  be  decomposed  into  a  global  mean  component  s,  a  periodic 
component  s,  and  a  random  motion  s'  (Reynolds  &  Hussain  1972,  Hussain  1983, 
1986).  The  total  variable  is  given  by 

s(x,  y,  z,  t)  =  s(x,  y)  +  s(x,  y,  t )  +  s'(x,  y ,  z,  t)  (85) 

A  precise  definition  of  the  periodic  and  random  components  is  as  follows.  The 
spanwise-averaged  lift  coefficient  as  a  function  of  time  Cx(f)  oscillates  quasi-periodically 
around  zero  due  to  vortex  shedding.  Each  period  of  oscillation  is  divided  into  Np  =  16 
segments.  The  first  segment  is  chosen  to  coincide  with  a  zero  lift  coefficient.  Over  a 
time  interval  0  <  t  <  T  which  includes  Nt  vortex  shedding  periods,  the  instants  cor¬ 
responding  to  a  particular  segment  k  are  denoted  by  for  n  =  1, 2, •  •  • , Nt-  Given 
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the  spanwise  average  <  s(x,y,t)  >  of  s(x,y,z,t),  the  spanwise-average  at  constant 
phase  k ,  also  called  a  phase-average,  is  computed  by  averaging  values  of  <  s(x,  y,  t)  > 
corresponding  to  identical  phases  in  the  time  series: 

1  nt 

<s(x,y)>k= -rr-Y,  <s{x,y,tkn)>  k  =  l,---,Np.  (86) 

n=l 

The  periodic  response  at  constant  phase  S*  is  obtained  from  the  Np  averages  at 
constant  phase  by  subtracting  the  spanwise-  and  time-averaged  component  s(x,y): 

s(x,y)k  =<  s(x,y)  >k  ~s(x,y)  (87) 

Once  the  periodic  component  is  known  at  each  phase,  the  random  fluctuation  s'  can 
be  found  from  (85).  By  construction,  the  time-averaged  periodic  component  s(x,y), 
as  well  as  the  time-  and  spanwise-averaged  random  fluctuations  s'(x,y )  are  zero.  The 
triple  decomposition  leads  to  the  definition  of  Reynolds  stresses  which  represent  the 
distinct  contributions  of  the  random  and  periodic  motions,  denoted  as  u-u'  and  utuj 
respectively,  and  related  by  the  equation: 

u'iu'j  =  U{Uj  —  TTiUj  —  UiUj.  (88) 
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Figure  1:  Real  part  of  the  inviscid  modified  wave  number 

Exact;  -  :  Fifth  order  upwind; - :  Seventh  order  upwind 

- :  4th  order  Pade;  •  •  •  :  4th  order  central 


Figure  2:  Imaginary  part  of  the  inviscid  modified  wave  number 
—  :  Fifth  order  upwind  biased; - :  Seventh  order  upwind  biased 


XiAxAt/c 


Figure  3:  Modified  wave  number  of  the  viscous  second  derivative  scheme 
-  :  Exact; - :  Sixth  order  central;  •  •  •  :  Second  order  central 


Figure  4:  Time  scheme  stability  diagram:  the  region  within  the  ellipsoid  is  unstable. 
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Chapter  3 

Validation  of  the  Numerical 
Method 


3.1  Introduction 

At  Reynolds  numbers  below  90,  the  flow  over  a  circular  cylinder  is  two-dimensional 
(Bloor  1964,  Williamson  1989,  1991a).  Its  mean  properties  in  the  near- wake  and 
on  the  cylinder  surface  have  been  extensively  documented  numerically  and  experi¬ 
mentally.  In  Section  3.2,  the  accuracy  of  the  numerical  method  is  examined  in  the 
computation  of  the  steady  near-wake  region  at  Reynolds  number  20.  The  compu¬ 
tational  domain  radius  required  to  minimize  the  effect  of  boundary  conditions  on 
near-wake  quantities  is  found  by  observation  of  the  pressure  distribution  in  the  wake 
to  be  300RC,  where  Rc  is  the  cylinder  radius.  Pressure  and  drag  coefficients  at  the 
cylinder  surface,  as  well  as  the  near-wake  velocity  distribution  are  found  to  compare 
well  with  experimental  and  other  computational  results. 

At  Reynolds  numbers  between  40  and  150,  two-dimensional  vortex  shedding  oc¬ 
curs  according  to  a  well  defined  Strouhal-Reynolds  number  relation  (equation  1). 
Section  3.3  presents  the  results  of  two-dimensional  simulations  at  Reynolds  numbers 
of  80  and  100.  The  computed  Strouhal  frequencies  are  in  good  agreement  with  the 
measurements  of  Williamson. 

The  computation  of  the  linear  stability  of  a  periodic  channel  flow  at  Reynolds 
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number  7, 500,  based  on  channel  half- width  and  centerline  velocity,  provides  a  three- 
dimensional  test  of  the  numerical  method.  Because  the  corresponding  perturbation 
equations  do  not  admit  three-dimensional  growing  modes,  the  perturbation  is  chosen 
as  the  slowest  decaying  Orr-Sommerfeld  mode.  The  error  in  the  perturbation-energy 
growth-rate,  computed  with  a  time  step  of  0.02/i/C/oo,  is  less  than  1%  of  its  exact 
value  on  a  mesh  containing  (129  x  32  x  32)  points  in  the  cross- channel,  streamwise 
and  spanwise  directions  respectively.  The  parameters  h  and  Uoo  above  refer  to  the 
channel  half-width  and  the  steady-state  streamwise  velocity  at  the  channel  centerline. 
The  validations  described  in  this  chapter  were  performed  with  the  fifth-order  accurate, 
upwind-biased  scheme. 


3.2  Steady  Cylinder  Flow 

At  a  Reynolds  number  of  20,  a  steady  recirculation  bubble  is  attached  to  the  cylinder 
surface.  The  salient  features  of  this  flow  are  the  bubble  length,  its  velocity  distri¬ 
bution,  the  drag  coefficient,  the  pressure  and  vorticity  distributions  on  the  cylinder 
surface  and  the  angle  of  separation  of  the  boundary  layer.  This  Reynolds  number  is 
chosen  because  its  properties  are  well  documented.  It  has  been  investigated  experi¬ 
mentally  by  Tritton  (1959),  Coutanceau  &  Bouard  (1977),  Thom  (1933)  and  Taneda 
(1956),  semi-analytically  by  Imai  (1951)  and  Nieuwstadt  &;  Keller  (1973),  and  nu¬ 
merically  by  Dennis  &  Chang  (1970),  Nieuwstadt  &:  Keller  (1973),  Takami  &;  Keller 
(1969)  and  Fornberg  (1980)  amongst  others. 

The  grid-converged  results  presented  below  were  obtained  on  a  computational 
domain  of  radius  300RC,  with  a  mesh  containing  (256  x  300)  points  in  the  radial  and 
azimuthal  directions  respectively  (Appendix  D).  Azimuthal  points  are  equispaced, 
while  the  radial  grid-spacing  increases  geometrically  with  a  stretching  factor  of  1.03. 

Table  (3)  summarizes  the  values  of  the  principal  flow  features  of  the  near  wake.  At 
Reynolds  number  20,  the  computed  drag  coefficient  (1.99)  is  within  3%  of  Tritton ’s 
experimental  value  of  2.05.  The  velocity  in  the  recirculating  region  is  very  small 
and  difficult  to  measure  experimentally.  Coutanceau  h  Bouard  (1977)  measured 
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Simulation 

Experiment 

Asymptotic 

Drag  Coefficient 

1.99 

2.05 

2.05 

Bubble  Length 

0.93  D 

0.93  D 

0.93  D 

Separation  Angle 

43.8  Deg. 

41.6  Deg. 

45.8  Deg. 

Minimum  Velocity  in  Bubble 

-0.031  Uco 

-0.040C/oo 

-OMOUoo 

Position  of  Velocity  Minimum 

0.42  D 

0.36  D 

0.43  D 

Table  3:  Steady-flow  topology  at  Reynolds  number  20 

a  minimum  velocity  in  the  bubble  of  -0.040f/<x,,  located  0.36D  downstream  of  the 
cylinder  surface.  The  computed  values  for  these  parameters  are  —  0.031f/oo  and  0.42D 
respectively.  They  are  within  3%  of  the  asymptotic  solution  of  Nieuwstadt  &  Keller 
(1973),  which  predicts  a  minimum  velocity  of  — 0.030f7c<>  at  0.43  D. 

The  separation  angle  of  the  boundary-layer  and  the  wall  vorticity  (figure  6)  com¬ 
pare  well  with  other  computational  and  experimental  results.  In  the  separated  flow 
region,  the  wall  vorticity  magnitude  is  smaller  than  Fornberg’s  prediction  (1980),  but 
in  good  agreement  with  the  results  of  Dennis  k  Chang  and  Nieuwstadt  k  Keller. 
The  wall  pressure  coefficient  (figure  5),  agrees  well  with  the  results  of  these  three 
computational  groups. 

Outside  the  recirculation  zone,  within  2.5  diameters  downstream  of  the  cylinder, 
experimental  and  asymptotic  streamwise  velocities  coincide  (figure  7),  and  differ  from 
the  computation  by  less  than  one  percent.  From  5  to  20  diameters  downstream 
(figure  8),  the  computed  streamwise  velocity  agrees  well  with  Imai’s  analytical  results 
(1951),  as  well  as  with  Nishioka  k  Sato’s  simulations  (1974),  but  not  with  the  results 
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of  Takami  &  Keller  (1969). 


3.3  Laminar  Vortex  Shedding 

Two-dimensional  simulations  in  the  regular  vortex-shedding  Reynolds  number  range 
establish  the  temporal  accuracy  of  the  numerical  method.  These  simulations,  at 
Reynolds  numbers  80  and  100,  further  test  the  robustness  of  the  far-field  boundary 
conditions  in  unsteady  flow.  Details  of  the  grid  construction  and  mesh  refinement 
studies  at  each  Reynolds  number  are  provided  in  Appendix  D. 

The  computed  Strouhal  numbers  are  0.152  and  0.164  at  Reynolds  numbers  80 
and  100  respectively  (figure  9).  These  agree  well  with  the  experimental  data  of 
Williamson  (1989),  given  by  equation  (1)  as  0.153  and  0.164  at  the  same  Reynolds 
numbers.  The  lift  coefficient  time  histories  at  both  Reynolds  numbers  (figures  93  and 
97,  Appendix  D)  demonstrate  that  the  flow  is  energetic  only  at  the  Strouhal  frequency 
in  the  near  wake.  Don  (1989)  suggested  that  numerical  boundary  conditions  could 
generate  spurious  modulating  frequencies  in  the  lift  response.  The  absence  of  such 
modes  in  the  present  simulations  indicates  that  inaccuracies  at  the  outer  domain 
boundaries  do  not  unduly  affect  the  solution  in  the  near  wake. 


3.4  Linear  Stability  of  a  Forced  Channel  Flow 

The  mathematical  formulation  of  the  linear  stability  of  a  periodic  channel  is  presented 
in  Appendix  E.  The  salient  feature  of  this  problem  is  the  forcing  introduced  to  drive 
the  flow  along  the  channel.  This  external  forcing  is  necessary  because  a  constant 
pressure  gradient,  corresponding  to  a  linear  pressure  distribution  along  the  channel, 
is  not  a  solution  to  the  compressible  boundary-value  problem,  and  cannot  be  used  as 
a  momentum  source.  In  the  forced  compressible  channel  formulation,  the  pressure  is 
constant  while  the  external  forcing  acts  as  a  source  of  momentum  in  the  streamwise 
direction. 

Appendix  E  presents  the  analytical  steady-state  channel  solution,  as  well  as  the 
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Mesh  size 

Time  step 

Decay  rate 

Error  in  decay  rate 

65  x  16  x  16 

1.00  x  lO"2 

-3.19  x  10"3 

-27.5% 

129  x  16  x  16 

0.97  x  10"2 

-4.92  x  10"3 

-11.8% 

129  x  32  x  32 

1.94  x  10"2 

-4.45  x  10"3 

1.0% 

Table  4:  Linear  stability:  computed  energy-decay  rates 

derivation  of  the  appropriate  small-disturbance  equations.  The  Mach  number  in  this 
stability  analysis  is  set  to  M =  0.1.  The  Reynolds  number  is  Re  =  7,500,  chosen 
to  match  that  in  the  work  of  Malik,  Zang  &  Hussaini  (1985). 

In  three  dimensions,  the  slowest  decaying  Orr-Sommerfeld  mode  is  chosen  as  the 
perturbation  of  the  steady-state  solution.  The  eigenvalue  of  this  mode  (equation  207) 
is  c  =  cr  +  ici  =  0.02962395  -  *0.00220154,  corresponding  to  a  perturbation-energy 
decay-rate  of  2c,  =  4.40  x  10-3. 

Three  simulations  are  performed  on  a  computational  box  of  size  (2 h  x  2vh  x 
2 irh)  where  h  is  the  channel  half-width.  The  three  meshes  contain  (65  x  16  x  16), 
(129  x  16  x  16)  and  (129  x  32  x  32)  points  in  the  vertical,  streamwise  and  spanwise 
directions  respectively.  Both  periodic  directions  have  equispaced  point  distributions, 
while  a  hyperbolic  tangent  with  a  stretching  factor  a  =  4.5  (equation  215)  is  used 
across  the  channel.  The  initial  perturbation  field  is  scaled  so  that  the  amplitude  of  the 
maximum  vertical  velocity  across  the  channel  is  10  4f/oo)  where  U<x,  is  the  steady-state 
channel-centerline  streamwise  velocity. 

The  equations  are  marched  forward  in  time  with  the  second-order  accurate  implicit 
scheme  described  in  Chapter  2,  using  a  velocity-based  Courant  number  of  0.025  for 
the  (65  x  16  x  16)  and  (129  x  16  x  16)  cases,  and  of  0.1  for  the  calculation  with 
the  finest  mesh.  All  simulations  are  run  for  2T0,  where  To  =  21.21  h/Uoo  is  the  time 
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necessary  for  the  perturbation  wave  to  propagate  through  the  streamwise  length  of  the 
computational  box.  The  computed  perturbation  energy  decay-rates  and  associated 
time-steps  on  each  grid  are  displayed  in  table  (4).  On  the  finest  mesh,  the  decay  rate 
is  predicted  within  1%  of  its  exact  value.  In  their  incompressible  two-dimensional 
computation  of  a  growing  perturbation,  Malik  et  a l.  used  a  Fourier  method  with 
4  collocation  points  in  the  streamwise  direction,  and  a  second-order  finite  difference 
scheme  on  a  Chebyshev  mesh  in  the  vertical  direction.  With  64  and  128  points  across 
the  channel,  the  errors  in  perturbation  energy  magnitude  after  two  periods  (t  =  2 T0) 
were  20%  and  4.7%  respectively. 


(p- Poo) /Ip 


Figure  5:  Re=20;  Wall  pressure  coefficient 

-  :  Present  results;  •  :  Dennis  &  Chang  (1970) 

x  :  Fornberg  (1980);  A  :  Nieuwstadt  &  Keller  (1973) 


Figure  6:  Re=20;  Wall  vorticity 

-  :  Present  results;  •  :  Dennis  &  Chang  (1970) 

x  :  Fornberg  (1980);  A  :  Nieuwstadt  &  Keller  (1973) 


Figure  7:  Re=20;  Rear  axis  streamwise  velocity 
Present  results;  •  :  Coutanceau  k  Bouard  (1977);  A  :  Nieuwstadt  k  Keller  (1973) 


Figure  8:  Re=20;  Rear  axis  streamwise  velocity  in  the  near  wake 

-  :  Present  results;  •  :  Imai  (1951);  x  :  Takami  k  Keller  (1969) 

o  :  Nishioka  k  Sato  (1974);  A  :  Nieuwstadt  k  Keller  (1973) 


Chapter  4 

Subgrid-Scale  Model  Performance 


4.1  Physical  and  Numerical  Parameters 

In  this  chapter,  results  from  three  large-eddy  simulations  of  the  wake  behind  a  circular 
cylinder  at  a  Reynolds  number  of  3, 900,  based  on  cylinder  diameter  and  free-stream 
velocity,  are  examined.  The  three  simulations  are  based  on  the  fifth-order  accurate, 
upwind-biased  scheme  for  the  convective  fluxes.  They  use  respectively  no  subgrid- 
scale  model,  the  fixed-coefficient  Smagorinsky  model  and  the  dynamic  model.  In 
the  reversed  flow  region,  mean  velocities,  as  well  as  mean  turbulent  and  periodic 
Reynolds  stresses  are  compared  to  experimental  data  obtained  using  a  Particle  Image 
Velocimetry  technique  (Lourenco  k  Shih  1993).  Downstream  from  the  recirculation 
zone,  within  the  first  10  diameters  of  the  wake,  mean  velocities  and  total  Reynolds 
stresses  available  from  a  hot-wire  experiment  (Ong  k  Wallace  1994)  are  used  for 
comparison.  Both  experiments  were  conducted  at  Re  ~  3, 900. 

The  grid  selected  for  the  simulations  is  described  in  section  C.3.4  of  Appendix  C. 
The  mean  velocities  and  total  Reynolds  shear  stress  obtained  on  this  mesh  using 
the  dynamic  subgrid-scale  model  are  grid  independent  in  the  first  ten  diameters  of 
the  wake  (Appendix  F).  The  time-step  used  in  all  simulations,  At  =  0.004Rc/t/oo, 
corresponds  to  a  free-stream  velocity  based  Courant  number  of  about  0.3.  Statistics 
were  compiled  over  approximately  6  vortex  shedding  cycles,  or  about  60RC/Uoo  time 
units. 
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The  remainder  of  this  chapter  consists  of  five  main  sections.  Section  4.2  examines 
the  issue  of  experimental  uncertainty  in  the  P.I.V.  and  hot-wire  measurements.  The 
impact  of  three-dimensionality  on  the  computed  flow  field  is  examined  in  section  4.3, 
where  in  addition  results  of  two-  and  three-dimensional  simulations  without  subgrid- 
scale  eddy  viscosity  models  are  compared.  Differences  between  the  three-dimensional 
simulations  are  discussed  in  section  4.4. 


4.2  Experimental  Parameters 

4.2.1  P.I.V.  Experimental  Uncertainty 


The  Particle  Image  Velocimetry  experimental  set-up  is  described  in  detail  by  Lourenco 
&  Krothapalli  (1988).  In  this  experiment,  the  cylinder  has  a  diameter  of  1.905  cen¬ 
timeters  and  is  39  centimeters  long.  The  free-stream  velocity  is  20.67  cm/s,  and  the 
Reynolds  number  based  on  diameter  and  free-stream  velocity  is  3, 900. 

The  cylinder  is  impulsively  accelerated  to  a  constant  velocity  in  a  towing  tank 
facility.  The  flow  in  the  wake  is  seeded  with  tracer  particles,  and  a  plane  in  the  wake 
is  illuminated  with  a  pulsed  laser  beam.  Photographs  of  the  pattern  projected  by  the 
tracers  are  taken  at  a  frequency  of  6 Hz,  from  which  velocity  information  is  extracted. 
The  data  in  the  recirculation  bubble  downstream  of  the  cylinder  is  obtained  from  93 
instantaneous  velocity-field  images  spanning  29  vortex  shedding  cycles.  These  data 
form  the  basis  for  the  comparisons  between  simulations  and  experiments  within  the 
first  3  diameters  of  the  wake. 

The  experimental  mean  velocity  and  Reynolds  stress  fields  provided  to  us  did 
not  display  the  expected  symmetries  and  anti-symmetries  about  the  wake  center- 
line.  Experimental  profiles  displayed  in  this  chapter  have  been  symmetrized  in  a 
post-processing  operation.  Denoting  mean  streamwise  and  vertical  velocities  by  u 
and  v  respectively,  and  symmetrized  fields  by  the  subscript  's',  the  post-processing 
operations  were: 


«.(£>  V)  = - o -  *?)  -  9 


(89) 
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(90) 


ufu's&v)  = 

vV3{(,ri)  = 

«V»(U)  = 


2 

>?)  +  t>VU, -??) 

2 

ttV.(£,  q)  -  «w(&  -n) 
2 


(91) 

(92) 


The  representative  symmetry  error  A /  in  a  symmetrized  quantity  /»(£,  i])  is  chosen 
as  the  maximum  error  at  each  location  £  normalized  by  the  maximum  of  /  at 


A/(0  = 


max,,  |/a(£,g) -/(£,??) | 
max,,  |/(£,»?)i 


(93) 


Symmetry  errors  in  Reynolds  stresses,  streamwise  and  vertical  velocities  are  displayed 
in  figures  (10)  and  (11).  Errors  in  the  vertical  mean  velocity  v  are  comparable  to  v 
itself  over  the  entire  domain  of  measurement.  Past  1  cylinder  diameter  downstream, 
symmetry  errors  in  the  streamwise  velocity  stand  at  approximately  5%  of  the  max¬ 
imum  local  velocity.  Global  Reynolds  shear  stresses  have  a  20%  error  for  x/ D  <  2. 
At  radial  locations  between  2.5 D  and  4 D,  the  error  is  around  30%. 


4.2.2  Hot-Wire  Experimental  Parameters 

The  experiments  were  conducted  in  the  windtunnel  of  the  Turbulence  Laboratory,  at 
the  University  of  Maryland.  The  test-section  of  the  windtunnel  has  a  rectangular  exit 
cross-section  of  1.2  x  0.7  meters.  The  freestream  velocity  was  4.2  meters  per  second, 
with  a  turbulent  intensity  of  0.7%.  The  Reynolds  number  was  3,900.  A  circular 
cylinder  of  1.43  centimeters  in  diameter  was  mounted  across  the  test-section  at  the 
tunnel  half-height  location,  approximately  7.3  meters  from  the  end  of  the  contraction. 

The  hot-wire  probes  were  operated  in  constant-temperature  mode  with  an  over¬ 
heat  ratio  of  1.35,  with  a  12-Channel  A.A.  Lab  Systems  Hot-wire  Anemometer  Sys¬ 
tem.  The  tunnel  velocity  was  monitored  using  a  pitot-static  probe  connected  to  a 
Barocel  Electronic  Manometer. 

Experimental  uncertainties,  provided  by  Ong  &  Wallace,  were  0.02  Uoo  on  mean 
velocities,  0.03  U on  Reynolds  stresses.  Symmetry  errors,  derived  by  post-processing 
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the  experimental  data  according  to  equations  89  through  92,  were  negligible  compared 
to  these  experimental  uncertainties. 


4.3  Impact  of  Three-Dimensionality  on  Mean 
Flow  Characteristics 

The  impact  on  mean  flow  quantities  of  three-dimensionality  is  examined  by  comparing 
the  results  of  two-  and  three-dimensional  simulations  which  use  no  subgrid-scale  eddy 
viscosity  model.  Both  calculations  were  performed  on  the  same  grid  in  the  plane  of 
mean  motion. 

At  the  cylinder  surface,  the  two-dimensional  simulation  yields  drag  and  rms  lift 
coefficients  of  1.42  and  1.74  respectively.  These  are  substantially  higher  than  the 
experimental  values  of  0.98  and  0.1  i  0.05  (Norberg  1987).  The  computed  three- 
dimensional  mean  drag  and  rms  lift  coefficients  are  0.96  and  0.07,  which  are  within 
experimental  uncertainty.  The  skin-friction  obtained  from  the  two-dimensional  sim¬ 
ulation  (1.33  x  10"2),  is  about  50%  higher  than  in  the  three-dimensional  calculation 
(0.87  x  10-2).  These  differences  in  magnitude  can  be  directly  observed  from  the 
time  history  of  the  lift  and  drag  coefficients,  displayed  in  figures  12  and  13.  The 
time-averaged  lift  coefficient  in  the  two-  and  three-dimensional  calculations  is  of  the 
order  of  10“3.  These  results  are  in  accordance  with  the  conclusions  of  other  re¬ 
searchers  (Braza  et  a i.  1986,  1990;  Tamura  et  a 1.  1990).  The  irregularity  of  the 
vortex  shedding  at  subcritical  Reynolds  number,  illustrated  by  the  modulation  of  the 
lift  coefficient  in  the  three-dimensional  calculation  (figure  12),  is  not  as  pronounced 
in  the  two-dimensional  simulation.  The  qualitative  behavior  of  the  lift  response  at 
Reynolds  number  3,900  is  similar  to  that  observed  by  Schewe  (1986,  figure  4. a,  page 
38)  at  Reynolds  number  2.64  x  105,  shortly  before  the  Reynolds  number  at  which  the 
drag  crisis  occurs.  Although  the  vortex  shedding  irregularity  in  the  three-dimensional 
flow  was  not  investigated  in  the  present  study,  some  factors  which  cannot  be  at  its 
origin  in  the  numerical  simulations  include:  free-stream  turbulence  in  the  oncoming 
flow,  the  effect  of  spanwise  end  plates  and  their  orientations,  the  influence  of  size  or 
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geometry  in  a  wind-tunnel,  and  an  interaction  between  the  vortex  shedding  process 
and  the  flow  in  the  wake  past  ten  diameters  downstream. 

The  most  outstanding  feature  of  the  two-dimensional  simulation  however  is  the 
absence  of  an  attached  recirculation  region  behind  the  cylinder  (figure  14).  The 
mean  centerline  streamwise  velocity  is  positive  at  all  stations  downstream  in  the 
two-dimensional  calculation  (figure  15).  The  wall  vorticity  distribution  (figure  16), 
indicates  that  separation  bubbles  at  the  cylinder  surface,  located  symmetrically  above 
and  below  the  wake  centerline  at  (108.1°  <  |0|  <  124.6°),  are  present  in  both  two  and 
three  dimensions.  Additionally  a  separation  bubble  is  attached  to  the  downstream 
face  of  the  cylinder  in  three-dimensional  flow.  Originating  at  6  =  85.3°  on  the  cylinder 
surface,  its  closure  point  lies  at  the  wake  centerline  at  L/D  =  1.56.  In  the  two- 
dimensional  simulation  the  pressure  drop  behind  the  cylinder  is  about  twice  as  large 
as  in  the  three-dimensional  case  (figure  17).  The  tangential  velocities  in  the  cylinder- 
surface  boundary  layer  are  correspondingly  greater  in  two  dimensions  (figure  18). 

These  differences  between  two-  and  three-dimensional  simulation  results  indicate 
that  three-dimensional  structures  strongly  influence  the  near-wake  at  Reynolds  num¬ 
ber  3, 900.  These  structures  consist  of  pairs  of  counter-rotating  streamwise  vortices 
(Hayakawa  &  Hussain  1989,  Bays-Muchmore  &  Ahmed  1993,  Mansy  &  aJ.  1994). 
Instantaneous  surfaces  of  constant  vorticity  norm  in  the  near-wake  (figure  19)  show 
that  these  structures  observed  experimentally  are  present  in  the  three-dimensional 
simulation.  The  cylinder  stands  at  the  bottom  of  the  figure,  the  flow  evolving  ver¬ 
tically  along  the  page.  The  detached  shear  layers  are  visible  at  the  top  and  bottom 
edges  of  the  cylinder.  The  7 tD  spanwise  extent  of  the  computational  domain  contains 
three  pairs  of  counter-rotating  streamwise  vortices.  This  is  consistent  with  the  obser¬ 
vation  of  Bays-Muchmore  &  Ahmed  that  there  is  about  one  pair  of  counter-rotating 
streamwise  vortices  per  cylinder  diameter  in  this  range  of  the  Reynolds  number.  In 
the  present  simulation,  these  vortices  are  present  between  the  two  primary  spanwise 
rollers  in  the  first  five  diameters  of  the  wake.  In  the  second  half  of  the  domain  (5  <  x/ 
D  <  10),  streamwise  structures  appear  to  be  shorter  and  more  diffused  than  in  the 
first  half.  This  point  is  further  examined  in  chapter  5. 
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4.4  Simulations  With  and  Without  Subgrid-Scale 
Models 

This  section  discusses  the  mean  flow  and  Reynolds  stresses  predicted  by  three  sim¬ 
ulations.  One  simulation  uses  no  subgrid-scale  eddy  viscosity  model,  a  second  one 
is  performed  with  the  standard  fixed-coefficient  Smagorinsky  model  (chapter  2),  the 
third  one  with  the  least-squares  version  of  the  dynamic  model.  Two  distinct  regions 
of  the  flow-field  are  examined  successively:  the  vortex  formation  region  (x/D  <  4), 
which  includes  the  cylinder  surface,  recirculation  bubble  and  recovery  zone,  and  the 
near  wake  from  4  to  10  diameters  downstream. 

In  the  following,  the  expression  ‘total  Reynolds  stress’  refers  to  the  sum  of  the 
periodic  (u,-  uj),  turbulent  (u't  «'•),  and  subgrid-scale  (r^)  components  of  the  Reynolds 
stresses. 


4.4.1  Cylinder  Surface 


The  main  features  of  the  mean  flow  at  the  cylinder  surface,  which  include  the  skin- 
friction,  back-pressure  and  drag  coefficients,  the  Strouhal  number  and  the  flow  sepa¬ 
ration  angles,  are  summarized  in  table  5. 

The  drag  coefficient  Cd  and  the  back-pressure  coefficient  Cpb  are  functions  of 
the  Reynolds  number  (Roshko  &  Fiszdon  1969,  Cardell  1993):  In  the  neighborhood 
of  Re  =  3,900,  the  drag  coefficient  increases,  while  the  back-pressure  coefficient 
decreases  with  increasing  Reynolds  number.  At  Reynolds  number  3,900,  their  values 
are  respectively  0.98  ±  0.05,  and  —0.90  ±  0.05. 


The  mean  drag  and  back-pressure  coefficients  are  determined  by  the  mean  veloc¬ 
ities,  viscous  and  Reynolds  stresses  in  the  wake:  letting  overbars  indicate  time  and 
spanwise  averaging,  the  mean  large-scale  streamwise  momentum  equation 
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integrated  on  the  rear  flow  axis  ( y  —  0) 


from  the  cylinder  surface  ( x/D  =  0.5)  to 
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infinity  leads  to 

d  /*oo 

Poo«L  +  Poo  -  P*/i?=o.5  +  -^J0  5{puv  +  pM-a12  +  Tu)dx  =  0  (95) 

Dividing  this  expression  by  ^u^/2  yields  the  back  pressure  coefficient  as 

Cp.  =  2  +  2^-  /  (p  u  w  +  p  u'v'  -  a u  +  r12)dx  (96) 

ay  J o.5 

where  all  quantities  are  referenced  to  the  freestream  values  poo?  u0 0  and  the  cylinder 
diameter.  A  momentum  balance  analysis  around  the  cylinder  relates  at  each  stream- 
wise  position  (x)  the  drag  coefficient  and  momentum  thickness  ( 0m )  to  turbulent 
intensities,  pressure  coefficient,  viscous  and  subgrid-scale  stresses: 

/-|-oo _  /*-{-oo  _  t  *4"Oo 

Cpdy  -2  J  undy  +  2  j  ^  (au  -  rn )dy  (97) 

Note  that  at  Reynolds  number  3, 900,  the  far- wake  momentum-thickness  Reynolds 
number  U^Bm/v  is  approximately  1,911,  since 

lim  em(x)  ~  l-CD  ~  0.49  (98) 

X-+CO  V  2 

Since  they  are  directly  related  to  spatially  integrated  mean  quantities  in  the  wake 
region,  Cd  and  Cpb  are  some  of  the  most  aggregate  quantities  of  the  flow.  As  such, 
they  are  relatively  insensitive  the  subgrid-scale  turbulence  model  used  in  the  large- 
eddy  simulations.  The  three  simulations  predict  these  coefficients  within  experimental 
uncertainty.  Similarly,  the  computed  surface  separation  angles  and  the  skin-friction 
coefficient  are  not  significantly  different  in  the  three  separate  calculations.  A  small 
difference  is  however  observed  in  the  predictions  of  the  Strouhal  number,  which  are 
5%  and  3%  below  the  experimental  value  of  0.215  in  the  dynamic  and  fixed-coefficient 
model  simulations  respectively. 

Overall,  the  observation  of  wall  statistics  reveals  insufficient  distinctions  between 
the  three  calculations  to  establish  a  performance  hierarchy.  We  turn  in  the  following 
section  to  a  more  detailed  examination  of  the  computed  mean  wake  statistics  for  this 
purpose. 
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4.4.2  Vortex  Formation  Region 

Recirculation  Zone 


The  analysis  of  differences  between  the  three  simulations  indicates  that  overall,  the 
dynamic  model  gives  the  most  accurate  prediction  of  the  mean  flow  in  the  recirculation 
region.  The  length  and  streamwise  velocity  distribution  in  the  mean  separation  bub¬ 
ble  are  close  to  experimental  observations  when  computed  with  the  dynamic  model 
(table  5).  In  the  calculations  without  model  or  with  the  fixed- coefficient  model,  the 
bubble  length  is  over-predicted  by  17%  and  29%  respectively.  The  prediction  with 
the  dynamic  model  (L/D  =  1.36)  is  within  two  percent  of  the  experimental  value. 

The  mean  streamwise  and  vertical  velocities  within  the  recirculation  bubble  are 
displayed  in  figures  20  and  21.  The  experimental  uncertainty  of  around  50%  (fig¬ 
ure  10)  on  the  bubble  mean  vertical  velocities  makes  comparisons  between  the  ver¬ 
tical  velocities  predicted  by  different  simulations  difficult.  The  streamwise  velocity 
profiles  across  the  wake  however  indicate  that  the  dynamic  model  simulation  offers 
more  accurate  predictions  than  simulations  with  no  model  or  with  the  fixed-coefficient 
Smagorinsky  model.  The  longer  recirculation  zone  in  the  latter  two  induces  a  down¬ 
stream  shift  of  the  entire  streamwise  velocity  distribution  on  the  rear  centerline  (fig¬ 
ure  15)  compared  to  the  distribution  obtained  with  the  dynamic  model.  This  shift  is 
apparent  in  the  profiles  shown  in  figure  20  at  stations  x/D  —  1.54  and  2,  where  the 
wake  is  wider  and  the  velocity  deficit  is  larger  in  the  fixed-coefficient  Smagorinsky 
model  simulation  than  in  the  computation  without  model. 

The  Reynolds  stresses  at  x/D  =  1.54  (figure  22)  illustrate  the  differences  between 
the  three  simulations  in  the  recirculation  zone.  The  dynamic  model  simulation  repro¬ 
duces  the  streamwise  stress  accurately  over  the  entire  width  of  the  wake  layer.  With¬ 
out  subgrid-scale  model,  the  peaks  of  the  streamwise  intensities  are  underpredicted 
by  18%.  In  the  fixed  coefficient  Smagorinsky  calculation,  the  streamwise  intensity  is 
damped  in  the  wake-layer  center  (—0.5  <  y/D  <  0.5):  its  peaks  lay  43%  below  their 
expected  value.  A  similar  pattern  appears  in  the  vertical  intensity  and  Reynolds  shear 
stress.  At  the  wake  center,  the  experimental  vertical  intensity  is  0.35  ±  0.04  U^.  The 
dynamic  model  prediction  (0.29  U^)  is  close  to  the  lower  experimental  bound,  while 
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the  simulations  without  model  (0.19  U £,)  and  with  the  fixed-coefficient  Smagorinsky 
model  (0.10  U £,)  underestimate  the  vertical  intensity.  Near  the  center  of  the  wake,  for 
-0.3  <y/D  <  0.3,  the  calculation  with  the  fixed-coefficient  Smagorinsky  model  gives 
a  Reynolds  shear  stress  of  the  wrong  sign.  All  simulations  predict  that  at  x/D  =  1.54, 
the  periodic  Reynolds  shear  stress  has  the  sign  opposite  of  the  turbulent  shear  stress 
near  the  wake  center.  The  fixed  coefficient  Smagorinsky  model  result  displays  a  large 
region  with  periodic  Reynolds  shear  stress  of  reversed  sign  (figure  23),  where  the 
maximum  amplitude  of  the  stress  (0.029  t/£>)  is  more  than  twice  that  obtained  with 
the  dynamic  model  (0.013  U £,). 

Recovery  Region 

In  the  recovery  region,  between  the  bubble  closure  point  and  x/D  =  4,  where  the 
flow  accelerates  despite  the  adverse  pressure  gradient  it  faces  (figure  15),  the  Reynolds 
stress  prediction  ability  of  the  three  simulations  is  qualitatively  identical  to  that  in 
the  recirculation  bubble.  The  periodic  Reynolds  stresses  at  x/D  =  2.5,  displayed 
in  figure  24,  are  underpredicted  by  simulations  without  model  and  with  the  fixed- 
coefficient  Smagorinsky  model.  The  dynamic  model  calculation  predicts  all  Reynolds 
stresses  within  the  experimental  uncertainty. 

Spanwise  intensities  are  not  available  experimentally,  but  the  predictions  of  the 
three  simulations  are  not  significantly  different  in  the  recirculation  region  except  near 
the  wake  center  (figure  22).  The  observable  differences  in  the  total  fluctuating  kinetic 
energy  (figure  25)  in  that  region  are  thus  largely  due  to  the  streamwise  and  vertical 
Reynolds  stresses.  The  dynamic  model  calculation  predicts  total  fluctuating  kinetic 
energy  levels  up  to  about  30%  and  90%  higher  than  those  obtained  without  model 
and  with  the  fixed- coefficient  Smagorinsky  model  respectively. 
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4.4.3  Near- Wake  Region 

Mean  Velocities  and  Total  Reynolds  Stresses 

In  the  near- wake,  between  4  and  10  D  downstream,  the  significant  differences  between 
the  three  simulations  in  the  predicted  mean  velocities  disappear  almost  entirely  (fig¬ 
ures  26  and  27).  Vertical  velocity  profiles  in  figure  27  are  displayed  at  x/D  =  3  and 
x/D  =  4,  the  only  two  locations  provided  by  Ong  k  Wallace  where  the  magnitude 
of  v/Uoo  is  substantially  larger  than  the  experimental  uncertainty  of  0.02Uoo.  The 
streamwise  velocity  at  x/D  =  4,  which  is  the  end  of  the  recovery  region,  is  more 
accurately  predicted  by  the  dynamic  model.  Aside  from  this  residual  effect  from  the 
recirculation  zone,  streamwise  and  vertical  velocities  are  within  experimental  uncer¬ 
tainty  at  all  downstream  stations. 

At  the  streamwise  locations  x/D  >  4,  the  experimental  uncertainty  on  the  total 
Reynolds  stresses  is  approximately  0.03£/£,.  The  three  simulations  yield  streamwise 
intensities  and  Reynolds  shear  stresses  within  the  experimental  uncertainty.  The 
vertical  intensity  is  significantly  underpredicted  by  the  fixed-coefficient  Smagorinsky 
model  calculation  at  x/D  =  4,  7  and  10  (figures  28  through  31).  These  figures  further 
indicate  that  the  differences  between  the  Reynolds  stresses  in  the  three  simulations 
diminish  with  downstream  distance. 

Time  and  spanwise  averaged  total  Reynolds  stress  contours  in  the  near-wake  are 
displayed  in  figure  32.  The  Reynolds  stress  distributions  at  Re  =  3, 900  are  qualita¬ 
tively  similar  to  those  at  Reynolds  number  1.4  x  105  (Cantwell  k  Coles  1983).  In  both 
cases  in  particular,  the  maximum  mean  streamwise  and  vertical  intensities  occur  near 
the  mean  bubble  closure  point,  away  from  and  on  the  wake  centerline  respectively. 
At  Re  =  3,900,  the  maximum  mean  streamwise  intensity  (0.1817^)  is  19%  lower  than 
at  Re  =  1.4  x  105.  This  difference  is  10%  in  the  case  of  the  maximum  mean  vertical 
intensity,  which  stands  at  0.39f7£,  at  Re  =  3, 900.  The  distribution  of  the  Reynolds 
shear  stress  inside  the  recirculation  bubble  cannot  be  compared  at  both  Reynolds 
numbers  because  it  was  not  provided  by  Cantwell  k  Coles.  However  the  maximum 
value  of  the  shear  stress  (0.1317^)  at  Re  =  1.4  x  10s  appears  to  be  only  4%  higher 
than  that  computed  at  Re  =  3, 900. 
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Model  Coefficients  and  Eddy  Viscosities 

The  time-averaged  dynamic  model  coefficient  C(x,y)  (equation  28)  scaled  on  the 
squared  Smagorinsky  constant  ( C 2  =  0.0652),  is  displayed  in  figures  33  and  34.  The 
ratio  of  the  fixed  Smagorinsky  coefficient  Cs  to  its  dynamic  counterpart  ijc(x,y)  is 
about  2  at  x/D  =  1,  2.5  at  x/D  =  5  and  3  at  x/D  =  10.  The  dynamic  model  coeffi¬ 
cient  within  the  turbulent  core  of  the  wake-layer  thus  generates  a  value  of  \f&{x~y) 
varying  between  0.13  and  0.2  in  the  first  10  diameters  of  the  wake. 

Mean  eddy  viscosities  obtained  from  the  fixed-coefficient  and  dynamic  model  cal¬ 
culations  are  displayed  in  figure  35  scaled  on  their  respective  maximum  values.  At 
Reynolds  number  3,900,  the  free  shear  layers  are  laminar  (Cardell  1993).  It  is  thus 
expected  that  the  eddy  viscosity  in  the  shear  layers  should  be  substantially  less  than 
further  downstream  in  the  wake,  a  behavior  reproduced  by  the  dynamic  model,  but 
not  by  the  fixed-coefficient  model  (figures  35,  36).  Another  difference  in  the  eddy 
viscosities  is  their  vertical  support:  it  is  bounded  by  the  wake  edge  at  all  downstream 
stations  in  the  dynamic  model  simulation,  but  not  in  the  fixed- coefficient  Smagorinsky 
model  case. 

Ten  diameters  downstream  of  the  cylinder,  the  dynamic  eddy  viscosity  is  about 
ten  times  larger  than  the  eddy  viscosity  generated  by  the  fixed-coefficient  model. 
The  resulting  total  viscosities  (u  +  vt)  are  different  by  a  factor  of  two  in  the  dynamic 
and  fixed-coefficient  model  simulations,  and  a  factor  of  three  in  the  calculations  with 
the  dynamic  model  and  without  model.  Despite  these  differences,  low-order  statis¬ 
tics  appear  independent  of  the  presence  or  absence  of  a  subgrid-scale  model  at  that 
location.  , 

Periodic  and  Random  Reynolds  Stresses 

The  convergence  of  the  results  from  all  three  simulations  toward  the  experimental 
data  in  the  near-wake  occurs  despite  the  decrease  with  downstream  distance  of  the 


contribution  of  the  periodic  motion  to  the  total  Reynolds  stresses  (figure  37),  mea¬ 
sured  at  each  downstream  location  by  the  ratio  (Matsumura  &;  Antonia,  1993): 

/-f  oo  _  /  r+oo  _  _ 
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The  periodic  component  of  the  motion,  a  spanwise-averaged  quantity,  is  two-dimensional: 
the  periodic  spanwise  intensities  are  nearly  zero  at  all  locations  inside  the  wake.  At  x / 

D  =  4,  the  periodic  contributions  to  the  streamwise,  vertical  and  shear  stresses  reach 
28%,  60%  and  26%  respectively.  Seven  diameters  downstream,  they  stand  at  19%, 
58%  and  26%.  Ten  diameters  downstream,  the  streamwise  and  vertical  periodic  com¬ 
ponents  continue  to  diminish,  with  contributions  at  13%  and  46%  respectively,  while 
the  periodic  Reynolds  shear  stress  maintains  an  approximately  constant  contribution 
level  of  27%:  within  the  first  10  diameters  of  the  wake. 

That  the  simulations  with  and  without  subgrid-scale  model  predict  similarly  ac¬ 
curate  Reynolds  stresses  at  x/D  =  10  is  thus  not  linked  to  the  dominance  of  periodic 
stresses  over  random  stresses.  It  is  rather  a  consequence  of  the  small  magnitude  of 
the  subgrid-scale  Reynolds  shear  stress  r12  relative  to  the  resolvable  shear  stress  u'v' 
in  the  wake,  an  effect  discussed  below. 

Subgrid-Scale  Shear  Stresses 

Profiles  of  the  subgrid-scale  shear  stress  magnitude  normalized  on  the  resolved  Reynolds 
shear  stress,  Ir^l/ln'  t>'|,  at  downstream  stations  x/D  =  3,  5  and  10,  are  shown  in  fig¬ 
ure  38  with  the  vertical  direction  y  scaled  on  the  wake  half-width  ( H ).  They  indicate 
that  the  relative  magnitude  of  the  subgrid-scale  shear  stress  within  the  turbulent  core 
of  the  wake  increases  with  downstream  distance.  This  behavior  is  expected  because 
of  the  coarsening  of  the  mesh  (see  Appendix  C)  and  of  the  conversion  of  large-scale 
energy  associated  with  the  organized,  periodic  motion  into  turbulent  kinetic  energy 
as  the  flow  evolves  downstream. 

The  magnitude  of  the  subgrid-scale  shear  stress  relative  to  its  resolvable  counter¬ 
part  is  small  however,  even  at  x/D  =  10.  The  overall  magnitude  of  subgrid-scale 
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stresses  scaled  on  resolvable  stresses,  measured  as 

r  Ir^l  dy /  r  \-^\  dy  (100) 

J- 00  /  J  —co 

is  displayed  in  figure  39.  The  subgrid-scale  shear  stress  norm  achieves  a  maximum  of 
5%  of  the  resolvable  shear  stress  at  x/D  ~  1.2.  From  4  to  10  diameters  downstream, 
the  same  norm  ranges  between  1%  and  2%. 

These  results  confirm  the  earlier  observation  that  the  impact  of  the  subgrid-scale 
model  on  the  mean  flow  is  small  in  the  computed  near-wake.  It  thus  appears  that  the 
mean  flow  is  principally  governed  by  the  largest  scale  motions  of  the  flow,  which  are 
accurately  resolved  on  the  grid  in  the  simulations.  This  conjecture  is  examined  in  the 
following  two  sections.  The  first  section  investigates  the  nature  of  the  small-scales 
present  in  the  near-wake  at  Re  =  3,900.  The  second  section  compares  experimental 
and  computed  one- dimensional  energy  spectra,  and  examines  the  impact  of  numerical 
dissipation  on  the  computed  fields. 

Small-Scale  Turbulence  at  Re  =  3, 900 

Fifty  diameters  downstream  of  the  cylinder  in  the  center  of  the  wake,  Uberoi  k  Frey- 
muth  (1969)  find  an  inertial  range  at  Reynolds  number  4,320  spanning  about  half 
a  wave-number  decade.  An  inertial  range  over  one  wave-number  decade  develops  at 
that  location  for  Reynolds  numbers  greater  than  approximately  15, 000.  In  direct  nu¬ 
merical  simulations  of  a  small-deficit,  time-developing  wake  at  a  momentum-thickness 
Reynolds  number  of  2, 000,  which  is  comparable  to  that  of  the  present  simulations, 
Moser  k  Rogers  (1994)  find  a  short  inertial  range,  spanning  about  half  a  wave-number 
decade;  Contours  of  instantaneous  spanwise  vorticity  visually  confirm  the  presence 
of  small-scale  turbulence  in  their  simulations. 

Experimental  velocity  spectra  at  Reynolds  number  3,900  provided  by  Ong  k 
Wallace  (1994)  indicate  that  a  fc_5/3  tangency  at  x/D  =  10  occurs  over  a  range 
of  similar  size  to  that  of  Moser  k  Rogers  (figures  40,  41).  In  the  £22 (kx)  spectra 
however,  the  vortex  shedding  frequency,  which  represents  a  large-scale  motion,  lies  at 
the  edge  of  the  j fc-5/3  range  at  locations  between  five  and  ten  diameters  downstream. 
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Furthermore,  the  high  wave-number  end  of  the  k~ 5/3  range  lays  at  approximately 
five  times  the  Strouhal  frequency  (figure  41),  a  number  lower  than  the  transition 
frequency  in  the  separated  free  shear  layers  at  Re  =  3, 900,  which  stands  at  six  times 
the  vortex  shedding  frequency  (see  Bloor  1964,  Norberg  1987,  Wei  &  Smith  1986). 

At  Reynolds  number  Re  =  5, 000,  visualizations  of  the  free  shear  layers  around 
the  cylinder  suggest  that  secondary  vortices,  formed  in  these  layers  and  entrained  into 
the  wake  at  about  two  diameters  downstream,  are  highly  organized  at  that  location 
(Ahmed  et  a 1.  1992).  At  Re  =  1,770  the  instantaneous  flow  pattern  between  30  and 
70  diameters  downstream  (Van  Dyke  1982,  p.  101),  appears  dominated  by  motions 
scaling  on  the  cylinder  diameter.  This  is  in  contrast  with  the  flow  at  Reynolds  number 
1.1  x  105  (Van  Dyke  1982,  p.130),  in  which  small-scales  are  clearly  present  in  the  near¬ 
wake  region. 

Based  on  the  evidence  above,  it  is  unclear  whether  small-scale  turbulence,  such 
as  that  observable  de  visu  in  the  simulation  of  Moser  &:  Rogers,  is  present  within  the 
first  ten  diameters  of  the  wake  at  Re  =  3, 900.  The  following  section  further  examines 
this  issue  by  presenting  pictures  of  the  computed  instantaneous  vorticity  in  the  wake, 
which  have  no  experimental  equivalents,  as  well  as  one-dimensional  energy  spectra. 

One-Dimensional  Energy  Spectra  and  Numerical  Dissipation 

Power  spectra  as  a  function  of  frequency  presented  in  this  chapter  were  calculated  with 
the  technique  used  by  Choi  &  Moin  (1990).  The  time  series  of  a  velocity  component 
g(x,  y,  z,  t )  is  defined  at  Nz  =  48  spanwise  points  on  the  interval  0  <  t  <  T.  The  series 
at  each  spanwise  location  is  considered  as  a  statistically  independent  realization.  The 
power  spectrum  of  the  function  g  is  then  calculated  as  the  average  of  the  spectra  of 
these  48  realizations.  Each  velocity  record  of  N  points  in  time  at  a  given  spanwise 
location  is  divided  into  m  time  intervals  of  length  Tm  =  2 T/(m  +  1)  with  a  50% 
overlap,  each  containing  M  =  2 N/(m  +  1)  points.  The  frequency  corresponding  to 
the  vortex  shedding  motion  is  wst  =  27 r/  where  /  is  the  dimensional  Strouhal  number 
(St  =  fD/U0 o).  The  maximum  resolvable  frequency  is  o>max  =  MTr/Tm,  while  the 
frequency  resolution  is  given  by  Auj  =  2tt /Tm. 
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The  power  spectra  presented  in  this  section  were  obtained  from  the  simulations 
using  the  dynamic  subgrid-scale  model  and  with  no  model.  The  data  was  collected 
at  88  points  over  a  time  interval  Wool D  =  30,  which  was  divided  into  m  -  7 
overlapping  segments  of  length  TJUoo/D  =  7.5  containing  22  points  each.  The  power 
spectra  were  thus  ensemble-averaged  over  7  x  48  =  336  velocity-trace  sets.  The 
maximum  resolvable  frequency  was  7.3  u>sti  with  a  resolution  of  0.67  ust • 

In  the  present  simulations,  the  grid  is  fine  enough  to  resolve  most  of  the  energy 
in  the  near-wake  (see  vertical  lines  in  figures  40,  41).  Small-scale  motions  present  in 
the  flow  are  thus  expected  to  be  apparent  in  the  simulation  data. 

Figures  42  and  43  respectively  show  the  instantaneous  vertical  vorticity  on  the 
symmetry  plane  y  =  0  for  5  <  x/D  <  10,  0  <  z  <  Lz,  and  the  streamwise  vorticity 
at  r/D  =  10  in  the  simulation  using  the  dynamic  model.  These  contours  reveal  that 
only  large-scale  features,  which  scale  on  the  cylinder  diameter,  are  present  in  the 
simulation. 

The  absence  of  small  scales  in  the  large-eddy  simulations  is  confirmed  by  the 
computed  one-dimensional  frequency  spectra  £ii(w)  shown  in  figure  44  compared  to 
the  experimental  results  of  Ong  k  Wallace  at  stations  x/D  =  5,  7  and  x/D  =  10. 
The  spectra  obtained  from  the  simulations  with  no  subgrid-scale  model  and  with  the 
dynamic  model  are  compared  in  figure  45  to  confirm  that  the  high-frequency  damping 
observed  in  the  large-eddy  simulation  is  not  attributable  to  the  filtering  operation. 
The  computed  spectra  discussed  below  are  those  obtained  from  the  dynamic  model 
simulation. 

The  computed  streamwise  spectra  agree  well  with  the  experiment  for  frequencies 
lower  than  approximately  2.5  u>st  at  x/D  =  5.  At  higher  frequencies,  the  computed 
energy  levels  fall  rapidly  below  their  experimental  values.  At  locations  further  down¬ 
stream,  the  maximum  well-resolved  frequency  decreases  as  the  grid  coarsens  in  the 
streamwise  and  azimuthal  directions.  Using  Taylor’s  hypothesis,  the  maximum  re¬ 
solvable  frequency  on  the  grid  is  10.5  u>st  at  five  diameters  downstream.  At  that 
location,  about  three-quarters  of  the  resolvable  frequency  range  are  thus  affected  by 
numerical  dissipation. 
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The  evolutions  of  the  computed  spanwise  wave-number  spectra  with  downstream 
distance  are  displayed  in  figure  46.  These  spectra  are  not  available  experimentally. 
The  energy  content  at  high  wave-numbers  increases  between  x/D  =  0.7  and  x/D  =  1 
in  the  three  spectra.  Past  the  bubble  closure  point  however,  the  energy  at  a  given 
wave-number  decreases  monotonically  with  downstream  distance  in  all  spanwise  spec¬ 
tra. 

That  the  influence  of  the  subgrid-scale  model  on  low-order  statistics  in  the  near¬ 
wake  is  small  is  thus  attributable  to  the  inability  of  the  numerical  simulation  to 
accurately  represent  frequencies  which  should  be  resolvable  on  the  grid.  The  observed 
resolved  frequency  range,  which  includes  only  a  quarter  of  the  frequencies  sustainable 
on  the  mesh,  makes  the  dynamic  subgrid-scale  model  operate  at  scales  which  contain 
negligible  amounts  of  energy. 

These  results  indicate  that  the  fifth-order  accurate  upwind  scheme  used  to  eval¬ 
uate  the  convective  derivatives  generates  enough  numerical  dissipation  to  affect  the 
turbulence  in  regions  where  the  mesh  is  fine  enough  to  resolve  the  mean  flow  velocity 
and  Reynolds  stresses.  Thus  although  this  scheme  appeared  to  be  a  good  candidate 
for  direct  simulations  (Rai  k  Moin  1991),  its  usefulness  for  large-eddy  simulations  is 
limited  by  the  fact  that  the  mesh  size  must  be  sufficiently  fine  everywhere  in  the  flow 
for  numerical  dissipation  to  be  negligible  compared  to  its  molecular  counterpart.  The 
consequent  grid-size  may  be  significantly  dictated  by  that  requirement,  as  the  case  of 
the  cylinder  near- wake  presented  above  indicates.  The  adequacy  of  upwind-biased  dif¬ 
ferentiation  schemes  to  large-eddy  simulations  being  in  question,  the  following  chapter 
examines  whether  the  impact  of  numerical  dissipation  on  the  computed  solution  can 
be  controlled  by  choosing  higher-order  accurate,  one-point  upwind  schemes  for  such 
calculations 

4.4.4  Simulation  Comparison  Summary 

Significant  differences  in  the  mean  velocity  and  total  Reynolds  stresses  between  the 
simulations  without  model,  with  a  fixed-coefficient  Smagorinsky  model  and  with  a 
dynamic  model  exist  in  the  vortex  formation  zone.  This  zone  consists  of  the  first  four 
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diameters  of  the  wake,  which  includes  the  cylinder  surface,  the  recirculation  bubble 
and  the  recovery  region.  The  dynamic  model  calculation  predicts  more  accurate 
mean  velocities  and  Reynolds  stresses  than  the  calculation  without  model.  The  fixed- 
coefficient  Smagorinsky  model  simulation  is  the  least  accurate  one. 

The  dynamic  model  gives  rise  to  an  eddy  viscosity  distribution  which  is  in  better 
agreement  with  our  expectations  based  on  the  flow  physics  than  the  fixed- coefficient 
model.  The  maximum  mean  eddy  viscosity  computed  with  the  dynamic  model  occurs 
in  the  wake  region  near  the  mean  bubble  closure  location.  On  the  contrary,  the  fixed- 
coefficient  Smagorinsky  model  generates  the  highest  level  of  eddy  viscosity  in  the 
separated  free  shear  layers,  which  are  laminar  at  Reynolds  number  3, 900. 

Downstream  of  the  formation  region,  between  four  and  ten  diameters  downstream, 
there  are  no  significant  differences  in  the  mean  velocities,  total  streamwise,  spanwise 
or  shear  Reynolds  stresses  obtained  in  the  three  simulations.  Differences  reside  prin¬ 
cipally  in  the  predictions  of  vertical  intensities,  which  are  more  accurate  with  the  dy¬ 
namic  model  than  without  model  or  with  the  fixed- coefficient  model.  However  these 
differences  between  the  three  simulations  diminish  with  downstream  distance.  The 
fifth-order  accurate,  upwind-biased  differencing  scheme  used  for  the  convective  terms 
appears  to  generate  enough  numerical  dissipation  to  overwhelm  the  contribution  of 
subgrid-scale  eddy- viscosity  models  in  coarse- mesh  regions  of  the  flow.  The  following 
chapter  examines  the  impact  of  numerical  dissipation  on  the  computed  flow-fields. 
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Simulation  Re  =  3, 900 

2-D  3-D  Fixed  Dynamic  Experiment 

no  model  no  model  coefficient  model 
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Figure  12:  Lift  coefficient  at  Re  =  3,900;  (a):  2D;  (b):  3D 

-  :  Total  lift;  •  •  •  :  Viscous  lift 

(Note  the  vertical  scale  difference) 


Figure  13:  Drag  coefficient  at  Re  =  3,900;  (a):  2D;  (b):  3D 
-  :  Total  drag;  •  •  •  :  Skin  -  friction 
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Figure  14:  Re  —  3, 900;  Mean  streamlines  i 

(a):  3D  no  model; 
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x/D 

Figure  15:  Re  =  3,900;  Centerline  streamwise  velocity  and  pressure  coefficient 

-  :  Dynamic  model; -  :  No  model;  •  •  •  :  Fixed  coefficient  model 

- :  2D  no  model;  x  :  Lourenco  &  Shih;  A  :  Ong  &  Wallace 
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Figure  16:  Re  =  3,900;  Wall  vorticity 

:  Dynamic  model; - :  No  model;  •  •  • :  Fixed  coefficient  model 

- :  2D  no  model 


0  45  90  135  180 

9 


Figure  17:  Re  =  3,900;  Wall  pressure  coefficient 

—  :  Dynamic  model; -  :  No  model;  •  •  •  :  Fixed  coefficient  model 

- :  2D  no  model;  x  :  Experiment  (Re  =  3,000;  Norberg 1987) 
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Figure  18:  Re  =  3,900;  Tangential  velocity  in  the  cylinder  surface  boundary  layer 
(a):  0  =  90°;  (6)  :  9  =  120°;  (c)  :  6  =  150° 

-  :  Dynamic  model; - :  No  model;  •  •  •  :  Fixed  coefficient  model 

- :  2D  no  model 
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Figure  19:  Re  =  3,900;  Instantaneous  constant  vorticity  contours  in  the  wake 
\wD/ud\  =  3;  Simulation  without  model 
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y/D  y/D 

Figure  22:  Re  =  3,900;  Total  Reynolds  stresses  at  x/D  =  1.54 
(a):  Streamwise;  (b):  Vertical;  (c):  Spanwise;  (d):  Shear 

-  :  Dynamic  model; - :  No  model;  •  •  •  :  Fixed  coefficient  model 

A  :  Lourenco  &  Shih 
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Figure  24:  Re  —  3,900;  Periodic  Reynolds  stresses  at  x/D  =  2.5 

-  :  Dynamic  model; -  :  No  model;  •  •  •  :  Fixed  coefficient  model 

A  :  Lourenco  &  Shih 


Figure  25:  Re  =  3,900;  Total  kinetic  energy  (u'k  u'k  +  Uk  uk  +  Tkk)/‘^U^0 
(a):  x/D  =  0.8;  (b):  x/D  =  0.9;  (c):  x/D  =  1.0;  (d):  x/D  =  2.0 
-  :  Dynamic  model; -  :  No  model;  •  •  •  :  Fixed  coefficient  model 
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Figure  26:  Re  =  3,900;  Streamwise  velocity  at  x/ D  —  4  (a);  7  (b);  10  (c) 

-  :  Dynamic  model; - :  No  model;  •  •  •  :  Fixed  coefficient  model 

A  :  Ong  &  Wallace 


y/D 

Figure  27:  Re  =  3,900;  Vertical  velocity  at  x/D  =  3  (a);  4  (b) 

■  :  Dynamic  model; - :  No  model;  •  •  •  :  Fixed  coefficient  model 

A  :  Ong  &  Wallace 


Figure  28: 
(c) 


Re  =  3,900;  Total  strearawise  Reynolds  stresses  at  x/D  =  4  (a);  7  (b);  10 

•  :  Dynamic  model; - :  No  model;  •  •  •  :  Fixed  coefficient  model 

A  :  Ong  &  Wallace 
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Figure  29: 
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Re  =  3,900;  Total  vertical  Reynolds  stresses  at  x/D  =  4  (a);  7  (b);  10  (c) 

-  :  Dynamic  model; - :  No  model;  •  •  •  :  Fixed  coefficient  model 

A  :  Ong  <fc  Wallace 
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Figure  30:  Re  =  3,900;  Total  spanwise  Reynolds  stresses  at  x/D  =  4  (a);  7  (b);  10 


:  Dynamic  model; - :  No  model;  •  *  •  :  Fixed  coefficient  model 

A  :  Ong  &;  Wallace 


y/D 


Figure  31:  Re  =  3,900;  Total  Reynolds  shear  stresses  at  x/D  =  4  (a);  7  (b);  10  (c) 

-  :  Dynamic  model; - :  No  model;  •  *  •  :  Fixed  coefficient  model 

A  :  Ong  h  Wallace 


Figure  32:  Re  =  3, 900;  Time  and  spanwise  averaged  total  Reynolds  stresses 
first  ten  diameters  of  the  wake 

(a):  Streamwise  stress,  contours  from  0  to  0.18,  increments  0.0075 

(b):  Vertical  stress,  contours  from  0  to  0.4,  increments  0.02 

(c):  Shear  stress,  contours  from  —0.12  to  0.12,  increments  0.01 


Figure  33:  Re  =  3,900;  Dynamic  model  coefficient  C(x,y)  scaled  on  the  squared 
Smagorinsky  constant  C]  at  xjD  —  0.5 
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Figure  34:  Re  =  3,900;  Dynamic  model  coefficient  C(x,  y)  scaled  on  the  squared 
Smagorinsky  constant  C] 


Figure  35:  Re  =  3,900;  Normalized  time  and  spanwise-averaged  eddy  viscosity  Vtj 
max(F<)  in  the  first  ten  diameters  of  the  wake 

(a):  Dynamic  model;  (b):  Fixed  coefficient  model 
Contours  from  0  to  1;  increments  0.05 


Figure  38:  Re  =  3,900;  Percentage  of  subgrid-scale  to  resolvable  Reynolds  shear 
stress 


x/D 

Figure  39:  Re  =  3, 900;  Percentage  of  vertically  integrated  subgrid-scale  to  resolvable 
Reynolds  stresses 

-  :  Streamwise;  •  •  •  :  Vertical; - :  Spanwise;  -  :  Shear 

Dynamic  model  simulation;  Re  =  3,900 


Figure  40:  Re  =  3,900;  Experimental  velocity  power  spectra  E\\  at  y  = 

k  Wallace,  1994);  •  •  •  :  x/D  =  5; -  :  x/D  =  7; - :  x/D  =  10; 

fc-s/3; - grid  cutoff  at  x/D  =  10 
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Figure  41:  Re  =  3,900;  Experimental  velocity  power  spectra  E22  at  y  = 

k  Wallace,  1994);  •  •  •  :  x/D  =  5; -  :  x/D  =  7; - :  x/D  =  10; 

£-5/3. - grid  cutoff  at  x/D  =  10 


x/D 

Figure  42:  Re  =  3, 900;  Instantaneous  vertical  vorticity  at  y  =  0,  5  <  x/D  <  10,  0  < 
z  <  Lz.  Dynamic  model  simulation. 

Contours  from  —5  to  6  by  0.3 
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Figure  43:  Re  =  3,900;  Instantaneous  streamwise  vorticity  at  r/D  =  10,  -5  <  y/ 
D  <  5,  0  <  z  <  Lz.  Dynamic  model  simulation. 

Contours  from  —5  to  3  by  0.3 


Figure  45:  Re  =  3,900;  One-dimensional  spectra  E\\  from  LES  and  DNS  at  (a) 
x/D  =  5;  (b)  :  xjD  —  7;  (c):  x/D  —  10 

- :  Dynamic  model;  •  •  •  :  No  model; 
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Figure  46:  Re  =  3,900;  One-dimensional  spectra  En,  E2 2  and  E33  at  y  =  0  from  the 
dynamic  model  simulation 

•  :  x/D  =  0.7; - :  x/D  =  1; - :  x/D  =  3;  :  x/D  =  5 

- :  x/D  =  7;  -  :  x/D  —  10 


Chapter  5 

Numerical  Dissipation  Impact  on 
the  Solution 


5.1  Motivation  and  Objective 

The  results  outlined  in  the  previous  chapter  indicated  that  the  numerical  dissipa¬ 
tion  generated  by  the  fifth-order  accurate,  one-point  upwind  biased  scheme  applied 
to  the  convective  fluxes  had  a  nefarious  impact  on  the  calculations.  This  scheme  is 
referred  to  below  as  scheme  A.  In  the  near-wake  region  between  four  and  ten  diam¬ 
eters  downstream  of  the  cylinder,  frequency  spectra  of  the  velocity  indicated  that 
a  substantial  range  of  frequencies  resolvable  on  the  grid  were  damped.  This  range 
widened  with  increasing  downstream  distance.  As  a  direct  consequence,  the  subgrid- 
scale  eddy-viscosity  models  were  probably  overwhelmed  by  numerical  dissipation  in 
the  near-wake. 

To  assess  the  impact  of  numerical  dissipation  on  the  computed  solution,  the  results 
obtained  with  the  fifth-order  accurate,  upwind-biased  scheme  must  be  compared  with 
those  generated  using  a  scheme  with  less  inherent  numerical  dissipation.  The  objective 
of  this  chapter  is  to  establish  such  a  comparison. 

Rai  &  Moin  (1991)  used  a  sixth-order  accurate  central  scheme  for  the  convective 
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fluxes  in  an  incompressible  turbulent  channel  flow  simulation.  The  solution  was  de- 
aliased  in  Fourier  space,  an  approach  not  well  suited  for  finite-difference  methods  in 
complex  geometry.  Further,  incompressible  Navier-Stokes  solvers  involve  no  divisions 
and  fewer  multiplications  of  primitive  variables  than  compressible  solvers.  The  sever¬ 
ity  of  the  aliasing  problem  may  thus  be  different  for  incompressible  and  compressible 
numerical  solutions. 

A  differencing  scheme  for  convective  fluxes  which  has  been  successfully  applied  in 
compressible  simulations  (Blaisdell  et  al.  1990,  Lee  et  a 1.  1992)  is  a  central  scheme 
in  which  derivatives  of  the  convective  fluxes  are  computed  as 
d  Id 

—pui(k\u  +  k2v  +  k3w )  =  +  k3v  +  k3w ) 

+  ^pUi-^-(kiU  + k2v  +  k3w) +  \(k1u  + k2v  +  k3w)-^-pUi  (101) 
2  ox  l  ox 

where  Arx,  k2  and  k3  are  geometric  coefficients  given  in  table  7  of  appendix  A.  This 
formulation  conserves  total  kinetic  energy  in  the  inviscid  limit  (Feiereisen  et  al.  1981), 
and  has  been  demonstrated  to  control  aliasing  errors  in  the  calculations  of  decaying 
homogenous  turbulence  (Blaisdell  et  al.  1990).  Applied  to  the  flow  over  a  circular 
cylinder  with  a  fourth-order  accurate  scheme,  this  method  was  found  to  be  unstable 
with  both  the  trapezoidal  and  the  second-order  backward  time-integration  schemes.  It 
was  discontinued  in  favor  of  a  strategy  to  systematically  reduce  the  level  of  numerical 
dissipation  generated  by  upwind-biased  schemes.  This  is  achieved  by  comparing  the 
velocity  power  spectra  obtained  from  simulations  with  increasing  order  of  accuracy 
which  use  one-point  upwind-biased  differencing  schemes  for  the  convective  fluxes. 

This  chapter  presents  the  results  of  such  a  computation,  which  uses  a  seventh- 
order  accurate,  upwind-biased  scheme  for  the  convective  fluxes.  This  is  referred  to  as 
scheme  B  below.  The  precise  differencing  stencil  is  given  by  equations  58  and  59  in 
chapter  2. 
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5.2  A  Numerical  Example 


The  solution  to  the  one-dimensional  wave  equation  provides  some  guidance  on  the 
improvements  which  may  be  expected  with  scheme  B.  The  governing  equation 

^  +  ^  =  0  0  <x<L,  t>  0  (102) 

at  ox 

with  the  associated  boundary  and  initial  conditions 

u(x,0)  =  exp(— a(x  —  xo)2)  u(0,t)  =  exp(— a(xo  +  t)2)  (103) 

are  advanced  in  time  using  the  second-order  accurate  implicit  scheme  used  in  cylinder 
simulations.  The  domain  length  is  chosen  as  L  =  126,  and  the  spatial  discretization 
is  uniform,  with  221  points  in  a  first  case  (1)  and  442  points  in  a  second  case  (2). 
The  examples  below  are  computed  with  the  parameters  x0  =  3  and  a  =  4,  which 
translate  into  a  Gaussian  of  half- width  1.07  spanned  by  four  and  eight  computational 
points  in  cases  (1)  and  (2)  respectively.  This  numerical  experiment  can  be  thought 
of  as  describing  the  effect  of  upwinding  on  a  coarsely  resolved  structure  traveling  on 
a  large-eddy  simulation  mesh.  The  solutions  using  schemes  A  and  B  are  compared  at 
the  time  when  the  traveling  Gaussian  is  centered  at  x  =  20.  The  exact  solution  is  a 
Gaussian  of  unit  amplitude,  unaltered  in  shape  from  its  initial  condition.  Figures  47 
and  48  display  the  profiles  of  the  computed  solutions  at  t  =  17.  Scheme  B  yields  an 
18%  improvement  over  scheme  A  in  the  amplitude  of  the  solution  with  four  points  per 
structure  (Gaussian).  With  a  doubling  of  the  resolution,  the  improvement  afforded  by 
the  seventh  order  scheme  diminishes  to  13%.  Simulations  with  both  schemes  signifi¬ 
cantly  damp  the  traveling  Gaussian  at  both  resolutions.  Both  solutions  also  exhibit 
the  effect  of  numerical  diffusion  by  widening  the  width  of  the  traveling  structure. 

This  one-dimensional  example  suggests  that  coarse  to  moderately  resolved  struc¬ 
tures  traveling  in  the  cylinder  wake  are  significantly  damped  by  the  numerical  dissi¬ 
pation  generated  by  upwinding,  although  the  application  of  scheme  B  should  result 
in  a  small  but  measurable  improvement  in  the  resolution  of  these  structures.  This 
conjecture  is  investigated  below. 
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5.3  Cylinder  Computations 

The  flow  at  Reynolds  number  3, 900  was  computed  on  the  same  (144  x  136  x  48)-point 
mesh  which  was  used  in  the  fifth-order  accurate  simulations  presented  in  the  previous 
chapter.  The  simulation  was  performed  with  the  dynamic  subgrid-scale  eddy-viscosity 
model. 

5.3.1  Instantaneous  Wake  Vorticity 

This  section  provides  a  qualitative  discussion  of  the  differences  observed  in  the  vortic¬ 
ity  fields  obtained  with  schemes  A  and  B.  In  all  figures,  the  vorticity  has  been  scaled 
on  the  local  mean  centerline  streamwise  velocity  deficit. 

The  vertical  vorticity  contours  at  r/D  =  5,  displayed  in  figure  49,  exhibits 
smaller  vortices  and  sharper  vorticity  gradients  obtained  with  the  higher-order  upwind 
scheme.  Ten  diameters  downstream  (figure  50),  similar  properties  are  discernable  in 
the  streamwise  vorticity.  In  the  calculation  with  scheme  A,  constant  streamwise  vor¬ 
ticity  contours  display  larger  structures  than  with  the  seventh-order  scheme. 

The  breakup  of  spanwise  vortical  structures  described  above  is  visible  in  the  span- 
wise  vorticity  contours  on  the  rear  plane  of  symmetry  ( y  =  0),  in  the  region  5  <  xj 
D  <  10  (figure  51).  The  simulation  with  scheme  A  displays  a  quasi-regular  pattern 
of  primary  rollers  spanning  the  entire  homogeneous  direction.  With  scheme  B,  these 
regular  rollers  have  been  altered  by  smaller  structures. 

Finally,  figures  52  and  53  display  the  instantaneous  streamwise  and  vertical  vor- 
ticities  respectively,  in  a  vertical  plane  at  z  =  0,  between  five  and  ten  diameters 
downstream.  Both  figures  exhibit  the  same  features  that  were  noted  in  the  discus¬ 
sion  above:  the  simulation  with  scheme  B  shows  smaller  and  more  numerous  vortical 
structures  in  the  wake,  separated  by  regions  of  sharper  vorticity  gradients. 

Differences  in  the  three-dimensional  structures  of  the  respective  wakes  are  fur¬ 
ther  illustrated  by  the  instantaneous  constant  vorticity  contours  (\uD/ud\  =  1,  3,  6) 
shown  in  figures  54  through  56  respectively,  where  ud  is  the  mean  centerline  velocity 
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deficit.  The  first  figure  offers  a  qualitatively  similar  view  of  both  simulations.  Stream- 
wise  structures  of  roughly  equal  size  appear  to  connect  the  spanwise  rollers.  However 
at  the  higher  levels  of  vorticity,  the  streamwise  ribs  are  thinner  in  the  higher-order 
accurate  simulation.  In  the  near-wake  (5  <  x/D  <  10),  streamwise  structures  appear 
as  slender,  coherent  ribs  extending  from  one  spanwise  vortex-core  to  the  next  with 
scheme  B.  On  the  contrary,  in  the  same  interval,  structures  computed  with  scheme 
A  are  shorter  ribs  which  diffuse  with  downstream  distance.  This  difference  is  most 
pronounced  for  \uD/u<i\  =  6,  the  corresponding  ribs  having  completely  disappeared 
from  the  near- wake  in  the  calculation  with  scheme  A.  This  indicates  that  the  ribs 
have  gained  in  strength  with  the  higher-order  accurate  scheme  and  lowered  level  of 
numerical  dissipation.  These  structures  are  instrumental  in  ejecting  low-speed  fluid 
from  the  wake-layer  and  entraining  irrotational  fluid  toward  the  turbulent  core  of  the 
wake.  Their  increased  strength  thus  corresponds  to  a  more  intense  mixing  process  in 
the  wake  layer. 

Instantaneous  contours  of  the  spanwise  vorticity  and  vertical  velocity  computed 
with  the  seventh-order  accurate  scheme  in  the  first  ten  diameters  of  the  wake  are 
displayed  in  figures  57  and  58  respectively.  The  first  two  Karmann  vortices  behind 
the  cylinder  form  a  pattern  similar  to  that  observed  at  laminar  regimes.  Past  five 
diameters  downstream,  the  breakup  of  the  spanwise  rollers  into  smaller  structures  is 
visible.  In  the  free  shear  layer  separating  from  the  bottom  side  of  the  cylinder,  the 
vertical  velocity  contours  outline  the  presence  of  secondary  vortices. 

5.3.2  Velocity  Power  Spectra 

The  observations  made  in  the  previous  section  on  the  presence  of  sharper  features  and 
smaller-scale  structures  in  the  seventh-order  accurate  simulation  are  quantified  in  this 
section  through  a  comparison  the  velocity  power  spectra  obtained  with  schemes  A 
and  B. 

At  locations  x/D  =  5,  7  and  x/D  =  10,  the  spectra  obtained  from  scheme  B 
(figure  59),  although  substantially  damped  by  numerical  dissipation,  display  energy 
levels  at  resolvable  frequencies  increased  by  up  to  about  eight  tenths  of  a  decade 
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relative  to  the  results  with  scheme  A.  Five  diameters  downstream  the  spectra  obtained 
from  the  seventh-order  accurate  calculation  are  damped  at  frequencies  lower  than 
that  of  the  start  of  the  k~ 5/3  range  displayed  in  the  experiment.  As  the  flow  evolves 
downstream  on  a  coarsening  mesh,  larger  scales  of  motion  are  expectedly  further 
damped.  To  accurately  resolve  these  spectra  at  xj D  —  5  past  the  start  of  the  inertial 
range  on  the  present  grid  would  thus  require  at  least  a  ninth-order  accurate,  one- 
point  upwind  scheme,  which  may  encounter  other  numerical  difficulties  associated 
with  Runge  phenomenon. 

Although  spanwise  wave-number  spectra  are  not  available  experimentally,  fig¬ 
ure  60  shows  that  the  gains  in  energy  at  high  wave-numbers  realized  with  scheme 
A  are  more  substantial  than  in  the  streamwise  frequency  spectra.  At  x/D  =  7  for 
instance,  the  energy  level  increased  by  up  to  two  decades  at  wave-numbers  between 
ten  and  fifteen. 

5.3.3  Dynamic  Subgrid-Scale  Model  Contribution 

The  magnitudes  of  the  subgrid-scale  intensities  and  shear  stress  relative  to  their  re¬ 
solvable  counterparts  are  shown  as  a  percentage  ratio  in  figure  61  over  the  region 
0.5  <  x/D  <  10.  In  the  near- wake  region,  the  streamwise  intensity  ratios  appear 
to  be  the  same  to  within  statistical  fluctuations,  whereas  it  is  higher  for  scheme  B 
beyond  x/D  ~  7.  In  the  same  zone,  the  vertical  intensity  ratio  is  relatively  constant 
with  both  schemes.  The  shear  stress  ratio  peaks  at  about  3.5%  in  the  simulation 
with  scheme  B,  which  is  100%  higher  than  the  peak  obtained  with  scheme  A  in  the 
near-wake  region.  At  locations  5  <  x/D  <  10,  the  magnitude  of  the  subgrid-scale 
shear  stress  relative  to  the  resolvable  shear  stress  is  larger  when  computed  with  the 
scheme  generating  the  least  numerical  dissipation. 

At  all  stations  downstream  of  the  cylinder,  the  magnitude  of  the  mean  eddy- 
viscosity  (figure  62)  is  larger  with  scheme  B.  A  maximum  gain  of  approximately  30% 
in  vt  occurs  at  x/D  =  5,  at  the  downstream  edge  of  the  vortex  formation  zone.  In  the 
near  wake,  differences  in  between  the  two  differencing  schemes  decrease,  so  that  at 
x/D  =  10  both  schemes  yield  comparable  distributions  of  eddy  viscosity. 


100 


5.3.4  Cylinder  Surface  and  Near- Wake  Mean  Flow 

The  statistics  of  the  flow  at  the  cylinder  surface  are  summarized  in  table  6.  The 
Strouhal  number  is  within  experimental  error  with  scheme  B  (St  =  0.217),  whereas  it 
layed  about  5%  below  Cardell’s  value  of  0.215  in  the  calculation  with  scheme  A.  The 
back-pressure  and  drag  coefficients  are  within  experimental  error  with  both  schemes, 
although  using  scheme  B  improves  their  prediction.  The  wall  vorticity  and  pressure 
distributions  predicted  by  both  methods  (figures  63,  64)  differ  slightly  on  the  wake¬ 
facing  side  of  the  cylinder,  where  scheme  B  yields  a  slightly  more  accurate  prediction 
for  the  back-pressure  coefficient  ( Cpb  —  —0.90)  than  scheme  A. 

In  the  vortex-formation  region,  the  centerline  velocity  is  lower  and  the  pressure 
coefficient  larger  with  the  higher-order  scheme  (table  6,  figure  65).  As  a  consequence, 
the  recirculation  bubble  is  longer  than  that  obtained  with  scheme  A.  In  this  region,  the 
turbulent  kinetic  energy  is  consistently  higher  with  scheme  B:  the  energy  present  in  the 
separated  free  shear  layers  is  several  times  larger  than  with  scheme  A,  indicating  that 
numerical  dissipation  has  been  significant  (figure  66).  At  x/D  =  1,  the  calculation 
with  scheme  A  predicts  a  higher  turbulent  kinetic  energy  in  the  bubble  core  than 
in  the  shear  layers.  With  scheme  B,  the  shear  layers  at  that  location  are  clearly 
delineated  and  contain  about  50%  more  energy  than  the  bubble  per  se. 

In  the  near-wake  region,  the  streamwise  velocity,  scaled  on  the  mean  centerline 
velocity  deficit  and  wake  half-width,  is  better  represented  by  method  B  across  most 
of  the  wake  layer  (figure  67).  Scheme  A  overpredicted  the  growth-rate  of  the  wake 
layer.  The  improvement  in  the  growth-rate  prediction  with  scheme  B  is  evident  in  the 
streamwise  and  vertical  intensities  (figures  68,  69)  which  peak  at  values  closer  to  the 
experimental  results.  The  Reynolds  shear  stresses  similarly  show  a  marked  improve¬ 
ment  with  scheme  B  (figure  70).  Ten  diameters  downstream,  the  shear  stress  profile 
predicted  with  scheme  B  is  in  good  agreement  with  the  experimental  measurements. 
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5.4  Comparison  Summary 


The  preceding  results  indicate  that  the  reduction  of  numerical  dissipation  levels  be¬ 
tween  the  fifth-  and  seventh-order  accurate  schemes  has  a  significant  impact  on  the 
computed  solution.  The  effect  of  switching  from  scheme  A  to  B  led  to  sharper  res¬ 
olution  of  the  transitioning  separated  free  shear  layers,  the  increased  energy  levels 
at  high  frequencies  and  wave-numbers,  as  well  as  the  thinning  and  strengthening  of 
the  rib  vortices  connecting  spanwise  rollers.  Improvements  in  the  solution  computed 
with  the  higher-order  scheme  are  equally  apparent  in  the  low-order  wake  statistics. 
Despite  these  ameliorations,  the  streamwise  one-dimensional  frequency  spectra  indi¬ 
cate  that  a  substantial  portion  of  the  resolvable  frequency  range  is  damped  even  with 
scheme  B.  Resolving  frequencies  up  to  the  inertial  range  past  the  vortex  formation 
zone  appears  to  require  either  a  higher-order  accurate  scheme  or  finer  resolution  of 
the  first  five  diameters  of  the  wake. 
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Figure  47:  Traveling  Gaussian,  4  points  per  structure 

- :  Initial  condition  ( t  =  0) 

7th  order  scheme  at  t  =  17; -  :  5th  order  scheme  at  t  =  17 


t 

i! 

j» 

i! 

i 

i  ! 
i  ! 

A 

i  ! 

i 

i  ! 
i  ! 

Figure  48:  Traveling  Gaussian,  8  points  per  structure 

- :  Initial  condition  (f  =  0) 

7th  order  scheme  at  t  =  17; -  :  5th  order  scheme  at  t  =  17 
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Figure  49:  Re  =  3, 900;  Instantaneous  vertical  vorticity  at  r/D  =  5 
Contours  from  —8  to  6  by  0.4 

(a):  7th  order  accurate  scheme;  (b):  3th  order  accurate  scheme 


x/D 

Figure  51:  Re  =  3,900;  Instantaneous  spanwise  vorticity  at  5  <  x/D  <10;  y  =  0 

Contours  from  —6  to  6  by  0.3 

(a):  1th  order  accurate  scheme;  (b):  5th  order  accurate  scheme 
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Figure  54:  Re  =  3,900;  Instantaneous  constant  vorticity  contours  \uD/ud\  =  1 
left:  1th  order-accurate  simulation;  right:  5th  order-accurate  simulation 
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Figure  55:  Re  =  3,900;  Instantaneous  constant  vorticity  contours  \ujD/ud\  =  3 
left:  1th  order-accurate  simulation;  right:  5th  order-accurate  simulation 
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Figure  58:  Re  =  3, 900;  Instantaneous  vertical  velocity  contours  in  the  first 
eters  of  the  vertical  plane  z  =  0  (seventh-order  scheme) 
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Figure  59:  Re  =  3,900;  One-dimensional  frequency  spectra  Eh(uj)  at  (a):  x/D  =  5; 
(b)  :  x/D  =  7;  (c):  x/D  =  10 

•  •  •  :  7tfeorder  scheme; - :  5t/lorder  scheme;  -  :  Ong  &  Wallace 
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Figure  60:  Re  =  3,900;  One-dimensional  wave-number  spectra  En(kz)  at  (a):  x/ 
D  =  5;  (b)  :  x/D  =  7;  (c):  x/D  =  10 

- :  7t/lorder  scheme;  -  :  5t/lorder  scheme 


Figure  61:  Re  =  3,900;  Subgrid-scale  Reynolds  stress  as  percentage  of  the  resolvable 
stress 

Streamwise  intensity  (a);  Vertical  intensity  (b);  Shear  stress  (c) 

- :  7f,lorder  scheme;  -  :  5i/lorder  scheme 
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Figure  62:  Re  =  3, 900;  Mean  eddy  viscosity  profiles 
(a):  x/D  —  1;  (b)  :  x/D  =  3;  (c):  x/D  =  4 
(d):  x/D  =  5;  (e)  :  x/D  =  7;  (f):  x/D  =  10 
- :  7thorder  scheme;  -  :  5t,lorder  scheme 


Table  6:  Cylinder  surface  and  bubble  region  comparisons 
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Figure  66:  Re  =  3, 900;  Formation  zone  total  kinetic  energy  (u'k  u'k  +  Uk  uk  +  Tkk)/2u2c 
(a):  x/D  =  0.8;  (b):  x/D  =  0.9;  (c):  x/D  =  1.0;  (d):  x/D  =  3.0 
-  :  5tfcorder  scheme; -  :  7'ftorder  scheme 
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Figure  67:  Re  =  3,900;  Streamwise  velocity  at  x/D  =  5  (a);  7  (b);  10  (c) 

-  :  5(,lorder  scheme; - :  7</lorder  scheme 

A  :  Ong  &  Wallace 
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Chapter  6 
Conclusions 


Three  large-eddy  simulations  of  the  flow  around  a  circular  cylinder  at  Reynolds  num¬ 
ber  3, 900  were  performed  using  the  dynamic  subgrid-scale  model,  a  fixed-coefficient 
Smagorinsky  model  and  no  subgrid-scale  model.  The  numerical  method  in  these 
three  simulations  was  based  on  a  fifth-order  accurate,  upwind-biased  finite-difference 
scheme  for  the  convective  terms,  a  sixth-order  accurate  central  scheme  for  viscous 
terms,  and  a  fully  implicit,  second-order  accurate  time  integration  technique.  Two 
simulations  using  the  dynamic  subgrid-scale  model  on  identical  grids  were  further 
compared;  one  was  based  on  the  fifth-order  accurate  scheme  described  above,  while 
the  other  one  used  a  seventh-order  accurate  upwind-biased  scheme  for  the  convec¬ 
tive  derivatives.  In  this  study  we  investigated  the  importance  of  three-dimensionality 
for  the  computed  solution,  the  performance  of  the  subgrid-scale  turbulence  models 
within  the  first  ten  diameters  of  the  wake,  and  the  impact  of  numerical  dissipation 
generated  by  high-order  accurate,  upwind-biased  numerical  schemes.  A  summary  of 
our  findings  is  presented  below. 

Three-Dimensionality  of  the  Near- Wake  at  Reynolds  number  3, 900 

The  near-wake  appears  to  be  strongly  three-dimensional.  In  a  two-dimensional 
simulation,  the  irregularity  of  the  vortex  shedding  present  in  the  three-dimensional 
solution  and  observed  experimentally  was  greatly  reduced.  The  three-dimensional 
near-wake  contains  counter-rotating  streamwise  vortices,  the  effect  of  which  cannot 
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be  reproduced  in  two  dimensions:  the  mean  two-dimensional  velocity  field  does  not 
reproduce  the  main  recirculation  region  behind  the  cylinder  which  is  present  in  the 
computed  three-dimensional  and  in  the  experimental  mean  fields. 

Performance  of  the  Subgrid-Scale  Turbulence  Models 

In  three-dimensional  calculations,  mean  wall  statistics  such  as  drag  and  pres¬ 
sure  coefficients,  wall  shear  stress  and  separation  angles  computed  with  the  three 
subgrid-scale  models  were  not  significantly  different  from  one  another.  In  the  vortex 
formation  region,  which  includes  first  four  diameters  of  the  wake,  the  calculation  with 
the  dynamic  model  predicted  more  accurate  mean  velocities  and  Reynolds  stresses 
than  those  without  model  or  with  the  fixed-coefficient  Smagorinsky  model.  Past  the 
vortex  formation  region,  differences  between  the  solutions  computed  from  the  three 
simulations  diminished  with  downstream  distance,  due  to  the  numerical  dissipation 
generated  by  the  one-point  upwinding  in  that  region. 

The  dynamic  model  predicted  a  mean  eddy  viscosity  distribution  in  better  qualita¬ 
tive  agreement  with  the  physics  of  the  flow  relative  to  the  fixed-coefficient  Smagorin¬ 
sky  model:  the  maximum  of  the  distribution  occured  in  the  turbulent  near-wake,  in 
the  vicinity  of  the  separated  bubble  closure  point.  With  the  fixed- coefficient  model, 
the  largest  mean  eddy  viscosity  was  generated  in  the  separated  free  shear  layers,  which 
are  laminar  at  Reynolds  number  3,900. 

Impact  of  Numerical  Dissipation  on  the  Computed  Solution 

The  seventh-order  accurate,  upwind-biased  scheme  generated  less  numerical  dissi¬ 
pation  in  the  wake  than  the  fifth-order  accurate  scheme.  As  a  result,  a  comparison  of 
the  simulations  using  both  schemes  revealed  that  the  magnitude  of  the  subgrid-scale 
shear  stresses  in  the  wake  increased  by  a  factor  of  about  70%,  and  that  the  energy 
content  of  the  velocity  fluctuations  was  larger  at  higher  frequencies  than  with  the  fifth- 
order  accurate  scheme.  Instantaneous  contours  of  the  vorticity  in  the  wake  provided 
further  evidence  as  to  the  presence  of  smaller-scale  structures  when  using  the  seventh- 
order  accurate  scheme.  Low-order  statistics  in  the  near-wake  showed  improvement 
with  the  higher-order  scheme,  although  streamwise  one-dimensional  frequency  spec¬ 
tra  were  damped  by  numerical  dissipation  over  a  significant  portion  of  the  resolvable 
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frequencies. 

The  fifth-order  accurate,  upwind  biased  scheme  had  shown  potential  for  direct 
numerical  simulations  of  flows  in  complex  geometries  (Rai  &;  Moin  1991).  Because 
of  the  numerical  dissipation  generated  by  the  scheme  on  coarse  meshes  however,  this 
method  was  found  to  be  unsuitable  for  the  large-eddy  simulation  of  the  wake  behind 
a  circular  cylinder.  In  the  near-wake  region,  where  the  mesh  was  coarse  but  yet  fine 
enough  to  predict  the  mean  velocity  and  Reynolds  shear  stresses  within  experimental 
uncertainty,  numerical  dissipation  annihilated  all  contributions  from  the  subgrid-scale 
eddy-viscosity  models. 

These  observations  suggest  that  the  adequacy  of  numerical  methods  for  large- 
eddy  simulations  cannot  be  extrapolated  from  conclusions  drawn  from  direct  numer¬ 
ical  simulations.  The  pros  and  cons  of  numerical  schemes  to  be  used  for  large-eddy 
simulations  should  be  investigated  using  large-eddy  simulations. 

Recommendations  for  Future  Research 

Within  the  context  of  large-eddy  simulation,  the  comparative  advantages  of  upwind- 
biased,  highly  accurate  methods  versus  central  schemes  need  to  be  assessed.  The 
former  control  aliasing  via  numerical  dissipation,  while  the  latter  can  be  formulated 
as  energy  conserving  schemes.  The  performances  of  these  methods  on  coarse,  highly 
stretched  grids  such  as  those  used  in  the  present  report  need  to  be  compared. 

Using  an  appropriate  numerical  scheme,  it  is  presumed  that  the  large-eddy  simu¬ 
lation  of  the  flow  around  a  circular  cylinder  at  Reynolds  number  1.4  x  105  (Cantwell  & 
Coles  1983)  would  reveal  substantial  performance  gaps  between  different  subgrid-scale 
models  in  the  near-wake,  because  of  the  presence  of  small-scales  in  that  region. 


Appendix  A 

Differencing  Schemes  and 
Generalized  Fluxes 


A.l  Differencing  Schemes  Near  Boundaries 

The  differencing  operators  applied  to  convective  and  viscous  fluxes  were  presented  in 
Chapter  2  (equations  56  through  61)  for  periodic  directions  and  wall-normal  directions 
away  from  boundaries.  Since  the  number  of  grid  points  available  to  differencing 
stencils  is  limited  by  the  presence  of  domain  boundaries  in  the  generalized  radial 
direction  £,  the  order  of  accuracy  of  the  corresponding  differencing  schemes  decreases 
near  the  cylinder  surface  and  far-field  boundaries.  The  mesh  spacing  in  computational 
space  is  set  to  Af  =  1.  Along  the  £  direction,  grid  points  are  ordered  from  1  to  N, 
the  first  point  being  located  at  the  cylinder  surface  by  convention.  The  following  two 
Sections  list  the  ^-derivative  operators  applied  to  the  convective  and  viscous  terms 
near  domain  boundaries. 

A.  1.1  Convective  Derivatives 

The  upwind  and  downwind  derivative  operators  are  denoted  by  and  super¬ 
scripts.  The  upwind  scheme  near  boundaries  is: 
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f  I,  -  -s* + £*  -  + f  L  -  ~F')+ *(A2) 

ttL-« Sf*-SF*+Sft+^+,(A*> 

TtL,  =  6SF"-3  -  +  ik^-1  +  jk*V  +  ,,<A3) 

- 1^- +  + 

A  second-order  accurate  central  scheme  is  used  at  boundaries.  Note  that  both  upwind 
and  downwind  fluxes  at  boundaries  are  evaluated  using  identical  central  schemes, 
which  is  equivalent  to  not  splitting  the  convective  fluxes  there.  The  downwind  scheme 


irL  - -JZF'-5KF2+iF°-ikF,+m3) 


dF^ 

N-2 


-sk^"-3  ~  5zFn-2 + iF"-'  _  ikFN + ^ 
%"L. =  2k {Fn  ~  f"-2) + ^ 


A.  1.2  Viscous  and  Heat  Flux  Derivatives 


Discretization  of  the  viscous  and  heat  fluxes  requires  first  and  second  derivative  op¬ 
erators.  Second-order  accuracy  is  retained  at  the  boundaries  in  the  corresponding 
numerical  operators.  First  derivative  schemes  applied  to  viscous  and  heat  flux  terms 
are  given  as  follows: 

I*  -  2k* + *<A2>-  f  L  -  >  - *>  ■ +  ”<A2) 

fl=i 
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dF 


|  =  _  F") +  3A^n_1  -  Fn~3^  +  ^(^4) 


dF 

dt 


1  f)F  1  O  Q 

^  =  i(F»  -  ft_)  +  «A*>,  ^  =  **-.  -  s*-.  +  +  ^(A2) 

The  second  derivative  schemes  are: 

S =  i(2Fi  -  + m  -  F<) + 1,(A2)’  l^L =  i(Fi  - 2F2 + Fs) + ,,(A2) 


d2F 


ae 


i^*  +  *>  +  E« +  *>-&*+«**> 


d2F 


a? 


W-2 


4 
3A 

_4 
3A 


d2F 


ati 


^^(/V-4  +  Fn)  +  J^(Fn~3  +  Fn- i)  -  2K^n~2  +  ^(A4) 

^  =  i( Fk.,  -  2Fh-i  +  fV)  +  -»(A2) 

ac2  at— i  A 

=  —(—Ftf- 3  +  4Fjv-2  ”  5Fv_i  +  2Fv) 


_  ,  .  <v 


A. 2  Inviscid  Fluxes  and  Flux- Vector  Splitting 


The  fluxes  F^,  Fv  and  Fz  in  directions  £,  tj  and  z  respectively  are  each  the  sum 
of  inviscid,  viscous  and  heat  fluxes,  Fe,  Fv  and  Fh  respectively  (see  Chapter  2). 
The  remainder  of  this  Appendix  details  the  expression  for  each  component  in  terms 
of  primitive  variables.  The  notation  in  this  section  is  identical  to  that  of  Beam  & 
Warming  (1978). 

Inviscid  fluxes  in  generalized  coordinates  are  given  by 


pu  \ 

pv  \ 

pu 2  +  p 

puv 

puv 

Xfj 

pu2  +  p 

puw 

pvw 

\u(e  +  p)  / 

\v(e  +  p)J 

(104) 
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pu2  +  p  puv 

F. !  =  -yt  puv  +  pu 2  +  p 


(105) 


u(e  +  p ) 


v(e  +  p) 


F%  =  J  pvw 
pw2  +  p 


(106) 


V  w(e  +  p)  / 

The  split  convective  fluxes  in  all  three  directions  can  be  written  as  one  expression 
by  introducing  geometric  coefficients  k\,  k2  and  k3.  The  general  form  of  a  split 
convective  flux  is 

/  2(7-l)Af  +  A±  +  A±  \ 

2(7  -  l)Af  u  +  A4  (u  +  cki)  +  A5  (u  -  cki) 

F*  -  2~  2(7  -  l)Af  v  +  Xf  (v  +  ck2)  +  A5  (u  -  ck2 )  (107) 

^  2(7-l)Afu>  +  X^(w  +  ck3) +  Xf(w  -  ck3) 

{ ^((u  +  cki)2  +  (f  +  ck2)2  +  (w  +  ck3 )2)  +  xi ) 


where 


(3  -  7)  fit  +  A5  )c: 

2(7  -  1) 


+  2p(7  —  l)Xf%i(k2w  —  k3v ) 


+  (7  -  l)Af  (u2  +  v2  +  w2)  +  -y((u  “  c^i)2  +  (v  “  CM2  +  (w  -  ck3 )2) 

and  where  the  superscript  ±  refers  to  the  construction  of  the  split  eigenvalues  A,  dis¬ 
cussed  below.  The  total  convective  flux  vector  is  the  sum  of  its  upwind  and  downwind 


components: 


F-t  =  F‘*  +  Pr,  f;  =  f;*  +  f;-,  f;  =  f*  +  P;-. 


(108) 


The  geometric  parameters  &,•  are  displayed  in  table  7  for  each  spatial  direction,  and 
ki  is  defined  by: 

b  =  ,  h - .  (109) 

'«  +  *?  +  *! 
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Direction 

kx 

k2 

k3 

£ 

Vv 

~xv 

0 

V 

XS 

0 

Z 

0 

0 

J 

Table  7:  Geometric  coefficients 


Direction 

Ai 

A4 

As 

* 

yvu  -  xvv 

V 

—y^u  +  X{V 

Ai  -Cyjxj  +  yl 

Z 

Jw 

Ai  -| -  Jc 

Ai  —  Jc 

Table  8:  Generalized  eigenvalues 


The  split  eigenvalues  A,  in  each  direction  are: 

_  A,  d:  |  At  | 


(110) 


where  the  convective  eigenvalues  A,  are  given  in  table  8.  The  split  eigenvalues  At 
defined  in  (110)  do  not  have  continuous  derivatives  in  A,-  at  sonic  and  stagnation 
points  (Van  Leer  1982,  Coirier  &  Van  Leer  1991).  This  singularity  leads  to  numerical 
oscillations  near  these  points  in  the  pressure  and  velocity  fields,  a  problem  which 
is  circumvented  by  fitting  A,-  with  a  parabola  in  the  vicinity  of  A i  =  0.  Let  Si  be 
the  threshold  value  of  |Ai|  below  which  the  split  eigenvalue  A,-  is  approximated  by  a 
parabola.  The  parabola  is  constructed  in  the  (A*,  A,)  plane  to  be  tangent  to  points 
(— 6,-,0)  and  (S{, Si),  and  in  the  (A,',  A,)  plane  to  be  tangent  to  points  (—Si,— Si)  and 
(<$;,  0).  The  modified  split  eigenvalues  are  then  given  by: 


A+  =  A,-  A,-  >  Si 

A?  1  1 

A+  =  — | — I-  -A,-  +  -Si  -Si  <  A{  <  Si 
4  Si  2  4 

A t  =  0  A ,  <  -Si 


(111) 

(112) 

(113) 
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and 


A,"  =0  \i>  Si 

-  A?  1  1 

K  =  -7r  +  o^i-76i  -Si<k<Si 

'  4  0;  2  4 

XT  =  A,  A,  <  -St 

In  a  single  formula,  these  expressions  reduce  to: 

^4  =  Ai±|Ai|  |A.|si.  Af  =  iA,±j(^+«i)  |Ai|<ft.  (117) 

In  practice,  we  choose  the  threshold  Si  to  be  the  largest  magnitude  of  two  eigen¬ 
values  A i  of  opposite  signs,  which  are  neighbors  along  a  grid  line.  Figure  (71)  dis¬ 
plays  the  modified  eigenvalues  X±  in  the  vicinity  of  a  singular  point  where  A  =  0. 
The  corner  formed  by  the  solid  lines  corresponds  to  the  singularity  generated  by 
flux-splitting.  The  eigenvalue  magnitudes  are  referenced  to  the  threshold  S ,  so  that 
parabolic  smoothing  is  applied  in  the  range  —1  <  |A|/<5  <  1.  A  quartic  smoothing 
function  was  also  tested,  and  produced  no  visible  differences  in  two-dimensional  test 
calculations  of  the  flow  over  a  circular  cylinder. 


(115) 

(116) 


A. 3  Iterative  Time-Scheme  Convective  Jacobians 

The  factored  iterative  equations  in  the  periodic  azimuthal  and  spanwise  directions 
(78  and  79)  involve  the  diagonalized  form  of  the  convective  flux: 

Ft  =  Pi\tPr1Q ■  (118) 


The  three  matrices  the  iterative  scheme  involves  are  Pv  x,  (Pz  1  Pv)  and  Pz,  which  are 
given  by  Pulliam  &  Chaussee  (1981)  as: 


-!/{( 1  -  (7  “  1  )$2c-2)  +  x^wp-1 
x4(l  -  (7  -  l)$2c-2)  +  y^wp-1 

P~l{Viv  ~  xtu) 

0((7  -  1)$V  -  c6) 

\  ^((7  -  1)$2^  +  c0) 


-yd 7  -  !)uc  2 

£f(7  —  l)uc~2 
xiP~x 

H-y^c  -  (7  - 1)«^) 
— /?(— t/^c  —  (7  —  1)«^) 
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(119) 


and 


-yd 7  -  Vvc  2 
x^(7  —  1  )vc~2 


-HP  1  -  vd'i  ~  i)™0  2  vdn  -  *)c  2  \ 
-ViP~l  +  x«(7  -  1  )wc~2  -xdl  ~  l)c-2 


where 


ViP'1 

0 

0 

P{ViC-(  7-  1W0 

1— H 

1 

1 

^(7-1) 

1 

o 

3 

1 

(7  -  i)vV>) 

/?0(7  —  l)u> 

^(7-1)  > 

0 

0  Vt 

Xt/y/2  -Xt/V2\ 

1 

0 

0  —X£ 

y,/V2  -yt/V 2 

p;1 

-Vi 

0 

0  0 

Y 

— 

xt/V% 

-yt/V 2  0 

0/2  0/2 

\  x^/y/2 

yt/V 2  0 

0/2  0/2  ) 

/  0  0 

1 

a 

a  ) 

0  - 

p 

u 

au 

au 

Pz  = 

P  0 

V 

av 

av 

0  0 

w 

a(w  +  c) 

a(w  —  c) 

\  pv  —pu 

$2  _|_  cwj 

a(*2  +  £-cw)) 

_ P_  . 

\/2  c 

$2  = 

\{u2  +  v2  +  w 2)  ; 

0  =  \Jx\  +  y| 

0  =  -j^u  +  x^v  ;  f3  = 


V2pc 
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A.4  Generalized  Viscous  and  Heat  Fluxes 


Viscous  and  heat  fluxes  F”’h  along  the  x  direction  are  decomposed  into  the  contri¬ 
butions  of  derivatives  of  the  primitive  variables  along  each  spatial  direction,  denoted 
as  F&,  Fvx  and  Fzx.  This  decomposition  is  presented  in  equations  (50)  through  (52). 
Expressions  for  each  of  these  fluxes  are  given  below.  In  each  vector,  the  symbol  E 
represents  the  dot-product  of  the  velocity  with  the  three  momentum  components. 
This  Section  completes  the  description  of  the  fluxes  which  appear  in  the  governing 
equations  of  motion  (55).  The  linearized  viscous  flux  Jacobians  in  the  left-hand  side 
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of  the  factored  equation  corresponding  to  the  wall-normal  direction  (80)  are  obtained 
from  the  viscous  fluxes  presented  below  by  straight-forward  differentiation.  Viscous 
and  heat  fluxes  are  functions  of  the  conservative  variable  vector  Q  and  its  spatial 
derivatives.  In  particular,  Fxx  =  FXX(Q,  Qx).  The  Jacobian  of  this  flux,  appearing  in 
the  left-hand  side  of  (80)  is  computed  at  time  level  n  +  1  from  the  Taylor  expansion: 

<K+1)  *  (^)“a<?"«  +  (^)”a Q»\  (122) 


The  remainder  of  this  Appendix  lists  the  expressions  for  the  viscous  and  heat  fluxes 
along  each  coordinate  direction. 


Viscous  fluxes  along  the  £  direction: 


*«  = 


JRe 


JRe 


(  0  \ 

(§»’  +  xV>ui  -  \xvVnvi 
(ixl  +  vDn  ~  \xvVr,n 
(xl  +  Vl)w( 
s  ) 

o 


-{\ynVi  +  xnxdun  +  (xnVt  ~  § Vvxt)vv 
~(iXVXi  +  yr,Vi)VV  +  (VriXi  -  lXr,ydUr> 
~ixixv  +  VvVt)wv 

E 

0 


f«  =  -rc 


-§  Vn*>z 
yvuz  -  xvvz 


Viscous  fluxes  along  the  //  direction: 

/  0 

-dVnVi  +  xnxdut  +  (Vvxi  ~  §**%)«€ 

~dxvx^  +  yr>ydvi  +  (xvVz  -  IwH 

~{xtxv  +  yvye)w(_ 

E 


Flr, 


1 

JRe 


(123) 
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(125) 
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Fv  — _ 
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Viscous  fluxes  along  the  z  direction: 


Fiz  ~  Re 
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Heat  fluxes  along  the  £  direction: 


(7  —  1  )PrReJ 
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(Xv+Vr,)-^ 


OT 
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Heat  fluxes  along  the  7  direction: 
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Heat  flux  along  the  z  direction: 


7  J  dT 
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Appendix  B 

Far-Field  and  Wall  Boundary 
Conditions 

B.l  Introduction 

This  Appendix  details  the  conditions  imposed  explicitly  and  implicitly  on  the  conser¬ 
vative  variables  at  the  wall  and  far-field  boundaries.  The  generalized  coordinates 
(£,  ?/,  z)  correspond  to  the  radial,  azimuthal  and  spanwise  directions  respectively. 
Lines  along  the  £  direction  extend  from  the  cylinder  surface  to  the  outer  bound¬ 
ary  of  the  computational  domain  (figure  84,  Appendix  F),  and  are  normal  to  the 
cylinder  boundary  surface  by  construction.  Wall  and  far-field  boundary  conditions 
are  required  at  the  £  boundaries.  The  azimuthal  ( rj )  and  spanwise  (z)  directions  are 
periodic,  and  the  factored  systems  in  these  directions  (equations  78  and  79)  are  solved 
with  periodic  boundary  conditions. 

B.2  Far-Field  Boundary 

Far-field  surfaces  are  partitioned  into  potential  and  wake  sub-regions,  different  sets 
of  boundary  conditions  being  appropriate  for  each.  In  potential  regions,  boundary 
conditions  are  based  on  Riemann  invariants  (Pulliam  1986b),  while  at  wake  bound¬ 
aries  the  primitive  variables  are  extrapolated  with  a  first-order  accurate  scheme  from 


their  values  inside  the  computational  domain.  The  locally-one- dimensional  Riemann 
invariants  R\  and  R2  are  defined  as 


Ri  =  Vn- 


R2  =  Vn  + 


2c 

7-1 


(135) 


where  c  is  the  local  sound  speed  and  Vn  the  velocity  normal  to  the  boundary  surface. 
Let  Vt  denote  the  velocity  tangential  to  the  surface.  According  to  our  axis  convention, 
in  which  the  7  direction  corresponds  to  the  azimuthal  direction  and  is  positive  as  in 
the  trigonometric  convention,  the  normal  and  tangential  components  of  the  velocity 


are  given  by  : 


Vn  = 


ynu  -  xvv 

\pn  +  Vl  ’ 


Xr)U  UrjV 

\]xl  +  y2v 


(136) 


B.2.1  Potential  Region  :  Explicit  Formulation 


In  the  following  derivations,  mesh  points  along  the  f  axis  have  indices  1  through  N , 
where  £1  and  lie  respectively  at  the  cylinder  surface  and  outer  boundary  of  the 
computational  domain. 

Inflow  boundary  conditions  are  chosen  as: 

—  —  (— )  ,  w  =  Woo j  Vt  =  Vtix„  R\  =  Rioo,  R2  =  R2N-1  (137) 

p 7  vp7/  00 

and  those  at  outflow  as 


P_ 

P1 


[  — )  ,  w  =  wn-  1,  Vt  =  K/v- 1)  R\  —  R\ooi  Ri  —  Rin-i- 

\pl  J  N- 1 


(138) 


Using  the  definition  of  the  conservative  variable  Q  (equation  44,  Chapter  2),  these 
relations  correspond  to  conditions  on  Q  at  the  boundaries  given  at  inflow  by: 


(139) 

(140) 

(141) 
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while  at  the  outflow: 


44  =  0 

qj  qj  +  ql  +  ql 

,s  7(1 -lji-'-1  2 


(142) 

(143) 


</l  =  — )7(R2N-1  ~  ^leo)2)  (144) 

\7piv-1  4  / 

42  =  ,f  ,  (*,vw.,  +  ijr,(Aloo  +  fijjv,,))  (145) 

\JXl  +  ylX  1  ' 

43  =  1  == - 2  (y^tN-l  ~  OO  +  R-2N-l))  (146) 
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(7-l)J^-i  '  2ft 
Conditions  at  infinity  are  found  from  the  potential  flow  around  a  fixed  circular  cylin¬ 
der: 
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B.2.2  Potential  Region  :  Implicit  Formulation 


The  time  scheme  being  written  in  delta-form  (see  equation  80),  implicit  boundary 
conditions  are  not  applied  to  the  conservative  variables,  but  to  their  incremental 
change  over  one  time-step  Aqt  =  A qf+1  =  q?+1  —  <?”.  Conservative  variables  at  far- 
field  boundaries  are  functions  of  the  invariant  il2,  the  pressure  p  and  the  tangential 
velocity  V*,  as  described  in  the  previous  Section,  so  that  their  increments  depend  on 
Af?2,  A p  and  AVt.  By  definition  these  increments  are: 
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Ap  =  (rzH('h+^ p±Aqt  _  P A?J  -  ®  A,3  -  — Ag*  +  A*) 

J  V  2ft2  ft  ft  ft  / 

AV$  = - /  !  f-— (xvq2  +  y„ft)Aft  +  *„Aft  +  yvA ft) . 

qi\/x2v  +  y2riK  qi  7 

Substituting  these  into  the  explicit  boundary  conditions  (139  through  148)  leads  to 
implicit  conditions  on  the  increments  A <7,.  The  resulting  implicit  inflow  boundary 


conditions  are: 
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while  at  the  outflow  boundary  implicit  conditions  are  given  by: 
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A <^5  =  $sA<7i  +  —  A<72  H - A93  H  A<74 
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where  the  parameters  $1  and  $5  are  defined  by 
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B.2.3  Wake  Region  :  Explicit  Formulation 

The  edge  of  the  wake  layer  at  the  outflow  boundary  is  defined  by  a  boundary-layer 
criterion 

u  <  0.99  uTO.  (160) 

Physically  realistic  boundary  conditions  within  the  wake  region  should  consist  of  4 
constraints  on  the  conservative  variables  and  take  the  viscous  stress  into  account 
(Poinsot  k  Lele  1989,  1992).  In  the  present  simulations,  outflow  conditions  are  sim¬ 
plified,  and  the  region  of  interest  in  the  computation  is  shielded  from  the  errors 
generated  by  inaccuracies  at  the  boundaries  by  a  long  region  in  which  the  mesh  is 
highly  stretched,  allowing  numerical  dissipation  to  damp  all  scales  of  motion.  At  the 
wake  boundary,  all  primitive  variables  p,u,v,w,p  are  extrapolated,  and  the  explicit 
conditions  are 


U  =  UN-1,  V  =  Vn-  1,  w  —  WN-l,  p  —  PN-l ,  P  =  PN-l- 


(161) 


B.2.4  Wake  Region  :  Implicit  Formulation 

The  implicit  boundary  conditions  corresponding  to  a  first-order  extrapolation  of  the 
primitive  variables  are: 

JN-i^qi  —  JA(<1i)n-i  (162) 
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A  p  =  Apiv-i 

(166) 

B.3  Solid  Walls 


B.3.1  Explicit  Formulation 


At  solid  surfaces,  no-slip  conditions  are  imposed,  =  93  =  94  =  0.  The  surface 
is  constrained  to  be  adiabatic  {d(q5/qi)/d£  =  0),  and  the  wall  pressure  derivative 
is  found  from  the  normal  momentum  equation.  The  normal  momentum  equation 
is  found  by  projecting  the  (  and  q  momentum  equations  onto  the  normal  direction 
vector  at  the  surface 

n  =  = — ={yvi-  xj),  (167) 

\Jyl  +  xn 


where  i  and  j  are  the  direction  vectors  in  the  x  and  y  directions  respectively.  The 
grid  is  normal  to  the  cylinder  surface  by  construction,  so  that  an  equation  governing 
the  wall  pressure  gradient  dp/d(  =  d(q5/J)/d£  is  necessary.  Since  the  time  scheme 
is  in  factored  form,  the  equation  defining  this  pressure  gradient  should  only  involve 
derivatives  in  the  (  direction.  The  appropriate  momentum  equation  is  obtained  by 
manipulation  of  both  the  normal  and  tangential  momentum  equations,  and  reduces 
to: 


dp  d2u  d2v  du  dv  d2u  d2v 

a'g( +  “2af  +  azW  +  °“s( +  “55{ +  +  ” 
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with  the  coefficients: 


<*i  =  x\  +  y2,  a2 


4  yv 


3  JRe 


144 


“6  =  +  Vixl  +  7x„x^„),  «7  =  ~3 jfej(8x*x*  +  x &n  +  7x^»)- 

Spatial  derivatives  of  the  pressure  and  temperature  at  the  wall  are  evaluated  with 
first-order  accurate  schemes. 


B.3.2  Implicit  Formulation 

The  no-slip  conditions  on  the  velocity  field  at  the  wall  hold  for  all  times,  so  that 

Aq-2  =  0,  Aq3  =  0,  A q4  -  0.  (169) 

The  implicit  boundary  condition  on  pressure  at  each  sub-iteration  neglects  the  change 
in  viscous  stresses  between  iterations;  the  implicit  adiabatic  wall  condition  is  obtained 
directly  from  its  explicit  counterpart: 

=  0,  A^r—  =  0.  (170) 

J  o(  qi 

The  boundary  condition  on  the  pressure  is  updated  after  each  time-step,  so  that 
neglecting  the  viscous  stresses  at  each  sub-iteration  does  not  affect  the  accuracy  of 
the  solution. 

B.4  A  Three-Dimensional  Example 

Figures  (72)  through  (76)  display  the  five  conservative  variables  at  a  given  instant  in 
time  in  a  large-eddy  simulation  of  the  turbulent  flow  around  a  circular  cylinder  using 
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the  dynamic  subgrid-scale  model.  The  Reynolds  number  is  3, 900  based  on  cylinder 
diameter  and  free-stream  velocity.  The  grid  contains  144  x  136  x  48  points  in  the 
radial,  azimuthal  and  spanwise  directions  respectively. 

The  figures  show  the  instantaneous  fields  in  the  entire  domain  (top  picture  on  each 
page),  with  the  outer  boundary  of  the  computational  domain  marked  by  a  limiting 
circle,  and  in  the  vicinity  of  the  cylinder  surface  (bottom  picture).  The  solution  near 
the  inflow  and  outflow  boundaries  does  not  exhibit  undesirable  numerical  oscillation 
in  any  of  the  conservative  variables.  From  the  contours  of  the  streamwise  and  vertical 
velocities,  the  numerical  dissipation  in  the  far  wake  of  the  large  scales  of  motion  is 
clearly  evident. 
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Figure  72:  Re  =  3, 900;  Instantaneous  density  ( p ) 
144  x  136  x  48  Dynamic  model  simulation 
top  :  whole  grid;  bottom  :  cylinder  vicinity 


Figure  73:  Re  =  3, 900;  Instantaneous  streamwise  velocity  (pu) 
144  x  136  x  48  Dynamic  model  simulation 
top  :  whole  grid;  bottom  :  cylinder  vicinity 


antaneous  vertical  velocity  (pv 
namic  model  simulation 
Dttom  :  cylinder  vicinity 


) 
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Figure  75:  Re  —  3, 900;  Instantaneous  spanwise  velocity  (pw) 
144  x  136  x  48  Dynamic  model  simulation 
top  :  whole  grid;  bottom  :  cylinder  vicinity 


Figure  76:  Re  =  3, 900;  Instantaneous  total  energy  (e) 
144  x  136  x  48  Dynamic  model  simulation 
top  :  whole  grid;  bottom  :  cylinder  vicinity 
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Appendix  C 

Cylinder  Grid  Generation 
Technique 


C.l  Introduction 

In  all  simulations,  the  computational  domain  is  an  annulus  of  outer  radius  RdRc  and 
spanwise  length  LZRC,  where  Rc  is  the  cylinder  radius,  spanned  by  (A^  x  Nn  x  Nz) 
points  in  the  radial,  azimuthal  and  spanwise  directions  respectively.  All  grids  share 
the  following  properties:  they  are  constrained  to  be  symmetric  in  the  upper  and 
lower  half  (£,  r])  planes;  generalized  radial  grid-lines  (r\  =  constant,  z  =  constant)  are 
normal  to  the  cylinder  surface  by  construction;  three-dimensional  grids  are  generated 
by  aligning  Nz  identical  two-dimensional  meshes  in  the  spanwise  direction;  spanwise 
points  are  equispaced;  azimuthal  points  are  distributed  in  the  (£,  7;)-plane  on  circles 
of  constant  radius  (£  =  constant). 

Given  these  characteristics,  the  grid  is  fully  described  by  the  location  of  each  of 
these  circles,  i.e.  the  distribution  of  radial  points  r(£),  and  the  layout  of  azimuthal 
points  along  each  circle  0{r]).  This  Appendix  describes  these  distributions,  which 
differ  for  laminar  and  turbulent  simulations.  The  following  two  sections  treat  these 
regimes  separately.  Figure  (77)  displays  the  nomenclature  used  throughout  this  report 
'  to  represent  the  features  of  the  near- wake  flow  structure.  The  locations  of  the  laminar 

surface  boundary  layer,  separated  shear  layers,  recirculation  and  potential  zones  are 
I 
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illustrated  using  contours  of  instantaneous  streamwise  velocity  extracted  from  the 
(116  x  136  x  48)  mesh  simulation  at  Reynolds  number  3, 900  presented  in  Appendix  F. 


C.2  Grids  for  Laminar  Simulations  at  Reynolds 
Numbers  20, 80  and  100 


C.2.1  Radial  Point  Distribution 


All  laminar  simulations  use  an  identical  family  of  grids,  described  in  this  Section.  In 
the  radial  direction,  the  mapping  from  physical  to  computational  space  is: 

r(£)  =  a  —  (1  +  a)e^ln4r  (171) 


where  a  is  a  constant  equal  to  1  +  (RD  -  1)/(1  -Sr(  *).  The  coefficient  sr  represents 
the  stretching  factor  of  the  grid.  Points  in  computational  space  are  equispaced  and 
given  by: 

6  =  *- 1  *  =  (172) 


The  mapping  r(£)  can  equivalently  be  written  as: 

r«)  =  1  +  (Ro  - 1)  * 

1  -  Sr 


(173) 


Substituting  the  expression  for  it  follows  that  rj+i— r,  =  A r,+i  =  (l+a)(l— sr)sj._1. 
The  discrete  grid  stretching  is  thus  defined  by  the  geometric  series  Ari+1  =  srAr, 
when  sr  is  constant.  Note  that  the  Jacobian  of  the  mapping  r(£)  is  a  smooth  function: 


dr  1  -Rd  £l 
^“l-^-lSr  r‘ 


(174) 


The  stretching  factor  sr  is  fixed  at  a  value  of  1.03  in  the  steady  flow  computations 
at  Re  =  20.  In  laminar  vortex  shedding  simulations,  the  coefficient  sT  is  itself  a 
geometrically  increasing  function  of  downstream  radial  distance: 

Sr(*)  =  si  +  («2  -  «i)— — *  =  1,  •••,-%  (175) 


153 


This  formulation  allows  for  slow  stretching  near  the  cylinder  boundary  characterized 
by  sj.  The  stretching  continuously  increases  in  the  wake  and  reaches  a  value  of  s2 
near  the  outflow  boundary.  Values  for  and  s2  are  given  in  tables  9,  10  and  11. 

C.2.2  Azimuthal  Point  Distribution 

Azimuthal  points  in  the  steady  simulations  are  equispaced: 

*(?)  =  2Xj i  (176) 

with 

Vj  =  i  -  1  j  =  (177) 

In  the  shedding  cases,  grid  points  are  selectively  distributed  in  two  separate  zones 
representing  the  potential  region  and  the  wake  layer.  The  grid  boundary,  or  envelope, 
of  the  wake  layer  in  the  plane  of  mean  motion  (x,y)  is  defined  analytically  by  a  line 
extending  from  the  cylinder  surface  to  the  domain  outer  boundary,  parametrized 
as  ( x(t),y(t ))  where  t  is  a  positive  continuous  parameter.  This  line  is  normal  to 
the  cylinder  surface  at  position  (cos  9q.  sin  9q),  and  reaches  the  outer  boundary  at 
(Rd  cos  $i,  Rd  sin^i),  where  Rd  is  the  computational  domain  radius,  90  and  9\  are 
free  parameters.  The  wake-layer  half  width  at  each  radial  position  is  y(t).  Within 
the  wake  layer,  a  given  number  A„wake  of  azimuthal  points  are  equispaced. 

Expressions  for  x(t)  and  y(t)  are  derived  by  considering  a  particle  traveling  in 
space  and  subjected  to  an  acceleration.  At  time  t  =  0,  its  position  and  velocity 
are  (cos  0O,  sin  0O)  and  (Uo,  U0  tan  90)  respectively.  The  location  of  this  particle  as  a 
function  of  time  is: 

(  x(t)  =  (Rd  cos  0\  —  cos  0O  +  Uo)t 2  —  Uot  +  cos0o  (178) 

\  y(t)  —  (Rd  sin  0i  -  sin  90  -f  Uo  tan  0O)*2  -  Uo  tan(0o)^  +  sin  90 

Although  this  envelope  could  be  used  to  define  the  wake  layer,  its  growth  downstream 
of  the  cylinder  ( y(t ))  is  only  proportional  to  t2.  However,  a  desirable  property  of 
the  grid  is  rapid  stretching  in  all  spatial  directions,  which  enhances  the  numerical 
dissipation  of  the  energetic  scales  of  motion,  downstream  of  the  part  of  the  wake  which 
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is  resolved  in  the  computation.  The  downstream  radial  distance  up  to  which  the  wake 
is  resolved  is  denoted  by  Rwake  in  the  following.  The  width  of  the  wake  envelope  can 
be  made  to  grow  exponentially  fast  provided  the  t 2  term  appearing  in  the  expression 
for  y(t)  is  replaced  by  (ebt  -  1  -  bt)/(eb  -  1  -  b)  where  b  is  a  free  parameter  which, 
together  with  Uo,  controls  the  width  and  growth-rate  of  the  envelope.  The  resulting 
expression  for  y(t)  is 

6^  —  1  —  bt 

y(t )  =  (Rd  sin  0\  -  sin  0O  +  U0  tan  0O)  b  _  ^  ^  ~  U0  tan(0o)t  +  sin  0O.  (179) 

The  grid  generation  parameters  used  at  Reynolds  numbers  80  and  100  are  dis¬ 
played  in  tables  (9)  and  (10).  In  simulations  at  Reynolds  number  80,  the  envelope 
equations  (178)  were  used  exclusively.  This  is  the  reason  for  which  no  values  of  the 
parameter  b  are  indicated  in  table  (9).  The  exponential  growth  of  the  wake  envelope 
is  a  feature  used  in  the  simulations  at  Reynolds  number  100. 

Figure  (78)  displays  the  wake  envelope  in  a  grid  constructed  with  parameters 
U0  =  —2,  b  =  10,  0o  =  85  and  6\  =  20.  The  same  envelope  in  the  near-wake  region 
is  shown  in  figure  (79).  The  wake  width  grows  as  the  square  root  of  radial  distance 
in  the  vicinity  of  the  cylinder,  and  exponentially  fast  in  the  downstream  wake. 


C.3  Grids  for  Turbulent  Simulations 

C.3.1  Azimuthal  Point  Distribution 

The  construction  of  computational  grids  for  the  turbulent  simulations  presented  in 
Chapter  4  and  Appendix  F  is  based  on  the  ‘wake  envelope’  described  in  the  previous 
section.  The  grid  is  equispaced  within  the  envelope.  The  one  exception  to  this 
approach  is  the  (174  x  128  x  48)  mesh  described  below.  In  that  grid,  the  azimuthal 
resolution  was  increased  by  50%  in  the  separated  shear  layer  regions  (40°  <  |0|  <  80°). 
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C.3.2  Radial  Point  Distribution 


In  the  radial  direction,  points  are  distributed  according  to  (171),  in  which  sr  is  a 
function  of  radial  distance,  defined  separately  in  three  regions.  Matching  at  the 
region  interfaces  ensures  that  the  mapping  from  physical  to  computational  space  has 
smooth  Jacobians.  These  three  regions  are  defined  by  5  parameters:  fixed  stretching 
ratios  at  the  cylinder  surface  and  outer  computational  domain  boundary,  (si)  and  (s2) 
respectively,  a  prescribed  maximum  grid  spacing  (Awake)  which  cannot  be  exceeded 
within  (/?wake)  radii  downstream,  and  the  grid-spacing  at  the  cylinder  surface  ( Arwau). 

In  the  first  zone,  the  grid  spacing  A r  grows  geometrically  from  its  wall  value 
Arwau  to  a  level  close  to  its  wake  value  Awake  at  a  rate  imposed  by  the  stretching 
Sj.  The  grid  then  becomes  equispaced  in  the  second  zone,  for  r  <  Rwake-  In  the 
outflow,  or  third  zone,  the  grid  spacing  increases  geometrically  from  Awake,  with  a 
stretching  ratio  which  increases  with  radial  distance  from  1  to  s2-  The  radius  of  the 
computational  domain  is  determined  by  the  total  number  of  radial  points  available. 

The  radial  increment  distribution  Ar(r)  is  displayed  in  figure  (80)  for  the  mesh 
containing  116  radial  points  with  the  corresponding  stretching  ratio  distribution: 

s,(ri)=iSf  <i8o> 

in  the  near-wake.  The  grid  spacing  is  Arwau  =  0.005i?c  at  the  wall  and  first  in¬ 
creases  geometrically  with  a  stretching  factor  sj  =  1.1.  In  the  resolved  wake  region 
extending  to  Rwake  =  10i?c,  the  increment  A r  is  close  to  the  prescribed  maximum 
value  of  Awake  =  0.25f?c  at  locations  past  2  cylinder  diameters  downstream.  Beyond 
the  resolved  region,  at  distances  r  >  i?Wake,  the  grid  stretches  geometrically  with  a 
stretching  factor  increasing  from  sx  =  1  to  a  maximum  of  s2  =  1.22.  In  this  example, 
the  unresolved  outflow  region,  extending  from  5  diameters  downstream  to  the  outer 
boundary  at  r  =  198.62/?c,  contains  26%  of  the  total  number  of  grid  points. 

C.3.3  Grid  Configuration  Summary 

This  section  describes  the  geometric  layout  of  each  mesh  used  for  turbulent  simu¬ 
lations.  The  relationship  between  grid  characteristics  and  physical  flow  features  is 
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examined  in  the  subsequent  Section. 

The  wake  envelope  and  near-wake  parameters  defining  the  grids  used  in  turbulent 
simulations  are  listed  in  tables  (11)  and  (12). 

The  mesh  refinement  study  presented  in  Appendix  F  involves  grids  containing 
(88  x  90  x  32),  (88  x  90  x  48),  (116  x  136  x  48)  and  (174  x  128  x  48)  points  in 
the  radial,  azimuthal  and  spanwise  directions  respectively.  The  spanwise  box  size  is 
Lz  =  2-kRc  in  all  cases.  The  corresponding  grid  configurations  in  the  near  wake  are 
displayed  in  figures  (81)  through  (83).  The  exponential  growth  of  the  wake  envelope  is 
illustrated  in  the  case  of  the  (116  x  136  x  48)  mesh,  shown  for  the  entire  computational 
(£,7 1)  plane  in  figure  (84). 

The  (116  x  136  x  48)  mesh  features  a  50%  uniform  increase  in  azimuthal  resolution 
and  a  100%  increase  in  the  radial  resolution  of  the  cylinder  boundary  layer  and 
resolved  near-wake  over  the  (88  x  90  x  48)  mesh.  The  mean  boundary-layer  velocity 
distribution  on  the  upstream  side  of  the  cylinder  is  found  to  be  identical  for  both 
simulations  on  meshes  with  90  and  136  azimuthal  points,  as  discussed  in  Appendix  F. 
The  finest  (174  x  128  x  48)  mesh  is  thus  constructed  with  the  same  azimuthal  point 
distribution  in  the  potential  region,  outside  the  ’wake  envelope’,  as  that  used  on  the 
(88  x  90  x  48)  point  mesh,  while  the  azimuthal  resolution  is  increased  by  50%  in  the 
separated  shear-layer  regions.  In  the  cylinder  boundary-layer  and  the  downstream 
recirculation  region,  the  radial  resolution  is  increased  by  100%.  The  resolved  wake 
region  is  extended  by  40%  to  7  diameters  downstream,  and  its  radial  resolution  is 
50  percent  finer  than  on  the  (116  x  136  x  48)  mesh.  Since  the  spatial  differencing 
scheme  is  globally  fifth-order  accurate,  mesh  refinements  by  50%  and  100%  in  a  spatial 
direction  represent  decreases  of  the  corresponding  leading  truncation  error  by  factors 
of  7.6  and  32  respectively. 

Section  F.4.1  of  Appendix  F  discusses  a  test  carried  out  to  evaluate  the  impact 
which  stretching  the  grid  radially  in  the  resolved  wake  region  has  on  the  mean  veloc¬ 
ities  and  eddy  viscosity.  The  (106  x  136  x  48)  grid  used  for  this  test  is  identical  to 
the  (116  x  136  x  48)  mesh  described  above,  except  that  between  1  and  5  diameters 
downstream  of  the  cylinder,  the  grid  stretches  geometrically  with  a  factor  of  1.02  in 
the  radial  direction.  Radial  increment  distributions  Ar(r)  for  the  unstretched  and 
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stretched  test-meshes  are  presented  in  figure  (85)  with  the  corresponding  stretching 
ratio  distributions  sr(r). 

The  final  grid  selected  for  large-eddy  simulations  contains  (144  x  136  x  48)  radial, 
azimuthal  and  spanwise  points.  This  grid  in  the  (£,77)  plane  is  displayed  in  figure 
(86). 

C.3.4  Physical  Flow  Features  and  Grid  Selection 

The  features  of  all  grids  used  in  turbulence  simulations  are  governed  by  four  principal 
length  scales:  the  cylinder  boundary-layer  thickness,  the  free  shear-layer  thickness, 
the  size  of  the  vortices  shed  from  the  cylinder  and  the  downstream  growth  of  the 
wake  layer. 

The  boundary-layer  thickness  on  the  upstream  side  of  the  cylinder  is  measured  by 
the  thickness  of  the  vorticity  layer  6U.  On  a  radial  line  (6  =  constant),  6W  is  defined 
as  the  distance  from  the  wall  at  which  the  vorticity  magnitude,  denoted  by  equals 
1%  of  its  maximum  value: 

fl2(r  —  6W)  =  0.01  max  fl2(r).  (181) 

1<t<Rd 

Because  the  surface  boundary-layer  separations  and  the  downstream  wake  have  a 
substantial  impact  on  the  pressure  distribution  on  the  cylinder  surface,  simplified 
momentum  integral  techniques  cannot  provide  accurate  estimates  of  the  vorticity 
thickness  distribution.  Approximate  values  of  6W  are  thus  obtained  from  the  most 
accurate  simulation,  performed  on  the  (174  x  136  x  48)  mesh.  Although  the  wake 
turbulence  is  not  fully  resolved  in  that  simulation,  the  mean  velocity  distribution  on 
the  upstream  face  of  the  cylinder  is  found  to  be  grid  independent.  Table  (13)  lists 
the  vorticity  thickness,  as  well  as  the  number  of  radial  points  in  each  grid  within 
the  vorticity  layer,  at  several  azimuthal  positions.  The  angle  6  in  that  table  has  its 
origin  at  the  downstream  mean  stagnation  point  ( x  =  Rc,y  =  0)  and  is  positive  in 
the  anti-clockwise  direction.  The  final  grid  selected  (144  x  136  x  48)  has  30  radial 
points  in  the  boundary  layer  at  the  top  of  the  cylinder  ( 6  =  90°),  and  about  14  points 
across  the  layer  at  the  forward  stagnation  point. 
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Cardell  (1993)  experimentally  documented  the  width  of  the  free  shear  layers  at 
Reynolds  numbers  4, 800  and  9, 500.  The  mean  width  6a  and  vertical  position  ys  of 
the  shear  layers  are  roughly  identical  close  to  the  cylinder  at  both  Reynolds  numbers. 
At  (x  =  0.5 D,ys  ~  0.6D),  6,  ~  0.15.D.  At  x  ~  2D  and  Re  =  4,800,  the  boundaries 
of  the  shear  layers  emerging  from  the  top  and  bottom  halves  of  the  cylinder  meet 
at  the  wake  center  and  6S  cz  0.8 D.  At  Re  =  9, 500  the  shear-layer  edges  meet  at 
x  ~  1.6D,  and  the  local  S3  is  also  about  0.8D.  These  data  lead  us  to  expect  that  at 
Re  =  3,900,  the  shear  layers  grow  from  about  6a  ~  0.15 D  at  (x  =  0.57),  ys  ~  0.67) ) 
to  6S  ~  0.87)  near  the  bubble  closure  point.  Interpolating  Cardell’s  data,  the  closure 
point  lay  at  x/D  =  1.83  ±  0.2  at  Re  =  3,900.  Table  (14)  lists  the  number  of  vertical 
grid  points  across  the  shear  layers  at  x/D  =  0.5  and  x/D  =  2.  The  (144  x  136  x  48) 
mesh  contains  9  and  16  points  across  the  shear  layers  at  these  respective  locations. 

The  number  of  grid  points  per  cylinder  diameter  in  all  three  directions  in  the  wake 
is  listed  in  table  (15).  In  the  (£,  rj)  plane  of  mean  motion,  the  numbers  shown  represent 
grid-point  densities  at  the  radial  limit  of  the  resolved  wake  region  (r  =  7?Wake)-  The 
azimuthal  resolution  at  that  location  is  the  coarsest  within  the  resolved  region,  since 
the  wake-grid  grows  with  the  wake  layer. 

At  Reynolds  number  3,900,  the  Kolmogorov  length  scale  (Ik)  at  locations  xj 
D  =  5,  7  and  10  are  displayed  in  table  16.  The  grid  spacing  in  Kolmogorov  units  ten 
diameters  downstream  for  the  final  grid  selected  is  27.2,  8.7  and  6.5  in  the  streamwise, 
vertical  and  spanwise  directions  respectively. 

The  wake-grid,  discussed  in  Section  C.2.2,  grows  downstream  with  the  square-root 
of  radial  distance.  The  adequacy  of  the  chosen  wake-envelope  growth  was  confirmed 
by  visual  inspection  of  the  flow  field.  Figure  (87)  displays  the  instantaneous  spanwise- 
averaged  streamwise,  vertical  and  spanwise  intensity  contours,  obtained  from  the 
(174  x  128  x  48)  simulation.  These  contours  are  superposed  onto  the  ((,y)  grid  and 
the  wake  envelope  is  clearly  seen  to  encompass  the  volume  of  active  turbulence. 
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Figure  77:  Diagram  of  near- wake  structure 
Snapshot  of  streamwise  velocity  contours  at  a  given  instant  in  time 
from  116  x  136  x  48  simulation  at  Re  =  3,900. 

Top  half  of  (£,  T})  grid  is  shown  with  skipped  grid-lines  in  radial  and  azimuthal  directions. 


Table  9:  Grid  construction  parameters  at  Re=  80 
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Grid 

Mesh  size  (r,  0) 

Si 

*2 

U0 

b 

00 

6i 

wake 

A  V  wall  /  Rc 

1 

129  x  81 

1 

40 

2 

129  x  121 

1 

80 

0.2078 

3 

129  x  141 

1.03 

l.i 

-5 

10 

85 

4 

200  x  141 

30 

100 

0.0250 

5 

256  x  141 

0.0048 

6 

300  x  171 

130 

0.0013 

Table  10:  Grid  construction  parameters  at  Re=  100 
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-50  0  50  100  150  2i 


Figure  78:  Wake  envelope  configuration 


r/Rc 

Figure  79:  Near-wake  envelope 
Construction  parameters:  Uo  =  —2,  b  =  10,  &o  =  85,  &\  =  20 
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Figure  80:  Near-wake  radial  spacing  (a)  and  stretching  ratio  (b)  distributions  for  the 
116  x  136  planar  mesh 
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Grid 

Mesh  size  (r,  0) 

U0 

b 

! 

00 

Sl 

s2 

Rd/Rc 

Purpose 

Number 

1,2 

88  x90 

-2 

10 

85 

20 

1.1 

1.22 

139.06 

Grid 

refinement 

study 

3 

116  x  136 

-2 

10 

85 

20 

1.1 

1.22 

198.62 

4 

174  x  128 

1.5 

10 

100 

20 

1.05 

1.20 

65.27 

Test  impact 
of  radial  grid 
stretching  in  wake 

106  x  136 

-2 

10 

85 

20 

1.1 

1.22 

176.07 

Final  grid  selected 

144  x  136 

-2 

10 

85 

20 

1.1 

1.22 

127.41 

Table  11:  Grid  parameters  for  wake  envelope  construction 
(Turbulent  simulations) 
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Mesh  size  (r,  0 ,  z ) 

N 

V  *7  wake 

Af  wall  /  Rc 

Rwake/  Rc 

Awake/  Rc 

88  x  90  x  32 
and 

88  x  90  x  48 

50 

0.01 

10 

0.5 

106  x  136  x  48 
and 

116  x  136  x  48 

75 

0.005 

10 

0.25 

174  x  128  x  48 

90 

0.0025 

14 

0.17 

144  x  136  x  48 

75 

0.0025 

21 

0.25 

Table  12:  Parameters  defining  the  resolved  wake  region 
(Turbulent  simulations) 
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Figure  81:  88  x  90  x  32  and  88  x  90  x  48  meshes  within  13D  of  the  cylinder 


Figure  82:  116  x  136  x  48  mesh  within  ISD  of  the  cylinder 
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the  cylinder 


Figure  84:  116  x  136  mesh  on  the  entire  domain 
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Table  13:  Radial  grid  points  in  surface  boundary  layer 


Location 

x/D 

Layer  thickness 

6./D 

Grid  size  (r,  6) 

88  x  90 

(116,106,144)  x  136 

174  x  128 

0.5 

0.15 

5 

9 

13 

2.0 

0.8 

10 

16 

21 

Table  14:  Vertical  grid  points  across  free  shear  layers 


Table  15:  Grid-point  density  in  the  near  wake 


x/D 

103  x  Ik/D 

Ax/Ik 

Av/Ik 

<1 

5 

8.81 

20.4 

6.4 

7.4 

7 

9.11 

24.1 

7.8 

7.2 

10 

10.1 

27.2 

8.7 

6.5 

Table  16:  Grid  spacing,  in  Kolmogorov  units,  in  the  final  grid  selected 
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Figure  87:  Instantaneous  spanwise-averaged  intensities  at  Re  =  3,900 
Illustration  of  the  wake  envelope  adequacy  using  the  (174  X  128  X  48)  simulation  data 
top:  Streamwise;  center:  Vertical;  bottom:  Spanwise 
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Appendix  D 

Laminar  Cylinder  Flow 
Validations 


D.l  Definitions 


This  section  defines  aerodynamic  coefficients  and  numerical  parameters  of  relevance 
in  cylinder  simulations.  The  pressure  and  viscous  components  of  the  drag  coefficient 
are  denoted  by  Cdp  and  Cdv  respectively,  and  those  of  the  lift  coefficient  as  Clp  and 
CLv .  The  total,  spanwise-averaged  lift  and  drag  coefficients  are  the  sum  of  these  two 
contributions: 


CL 


Lift 


\Po 0U12RCL; 


=  CLp  +  Civ,  Cd  = 


Drag 


\PooUl2RcL, 


=  CDp  +  CDv  (182) 


with 


1  fLz  [2ir  /A  \ 

CLv  =  D  ,,,  r  /  /  (  -(«4  cos  e  +  V£  sin  9)  sin  6  +  DZo  cos  6 )  dOdz 

ReM^Lz  Jo  Jo  \ 3  / 

1  rLz  f2* 

CDp  =  -  T7-2  -r  /  /  pcosO  d9dz 

MIqLz  Jo  Jo 
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Cdv  =  p  ,  /  /  ( cos  0  +  Vi  sin  6)  cos  9  —  QZq  sin  o) dOdz 

ReM^Lz  Jo  Jo  \S  / 

where  flZo  =  is  the  vorticity  at  the  wall  and  Lz  is  the  spanwise  box  size. 

The  Mach  number  M ^  appears  in  these  definitions  because  the  reference  velocity  was 

chosen  as  the  sound  speed  (Chapter  2). 

The  Courant  numbers  used  in  this  report  are  based  on  free-stream  velocity  ( CF L(Uoo)) 
or  sound  speed  (CFL(c0 0))  and  are  defined  as: 

lul  +  c 


CFLfcoo)  =  max 


(x,y,z)  L  Ax 


CFL(Uoo)  =  max 

(x,y,z) 


+ 


ju|_ 

Ax 


t,|  +  c+H+clA( 

Ay  A  z  J 

(183) 

M  Mi .  , 

^y  +  Ai\At 

(184) 

where  the  ‘max’  operator  extracts  the  maximum  value  of  its  argument  in  the  entire 
computational  domain. 


D.2  Steady  Flow  at  Re  =  20 

At  Reynolds  number  20,  the  flow  around  a  circular  cylinder  is  steady,  and  a  pair  of 
vortices  is  attached  to  the  downstream  side  of  the  cylinder.  This  Section  presents 
a  study  of  the  impact  of  domain  size  on  the  steady  solution.  It  is  shown  that  a 
small  computational  domain  does  not  affect  the  velocity  field  in  the  vicinity  of  the 
cylinder,  but  that  the  pressure  in  the  wake  region  is  sensitive  to  the  position  of  the 
outflow  computational  boundary.  An  appropriate  domain  radius  for  which  the  near¬ 
wake  pressure  and  velocity  are  not  affected  by  the  boundary  conditions  is  found  to 
be  about  300  cylinder  radii.  The  steady  solution  on  that  domain  is  compared  to 
experimental  results  in  Chapter  3. 

D.2.1  Effect  of  Computational  Domain  Size 

The  computational  domain  is  a  circle  of  radius  Rd-  In  all  simulations,  the  radial  grid 
distribution  stretches  geometrically  from  the  cylinder  surface  to  the  outer  boundary 
with  a  stretching  factor  fixed  at  1.03,  while  azimuthal  grid  points  are  equispaced. 
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Asymptotically  far  downstream  of  the  cylinder  the  velocity  defect  on  the  wake 
centerline  decreases  as 

« ~  -4-  •  (i85> 

Vr 

The  effect  of  the  body  is  thus  felt  at  large  distances  downstream,  so  that  extending 
the  computational  domain  to  distances  where  the  influence  of  the  cylinder  is  negligi¬ 
ble  is  not  practical.  Since  one  cannot  impose  exact  outflow  boundary  conditions  in 
the  wake  region,  the  mismatch  between  flow  physics  and  numerical  constraints  cre¬ 
ates  spurious  pressure  waves.  However  the  mesh  stretching  near  the  outer  boundary 
increases  numerical  dissipation,  dampening  these  perturbations.  The  influence  of  the 
boundary  conditions  on  the  near  wake  is  demonstrated  on  three  different  domains 
of  size  Rd  =  60i?c,  120 Rc  and  300i?c.  The  solutions  are  grid  independent  on  each 
domain  studied.  For  instance  on  the  larger  domain,  simulations  are  performed  with 
mesh  sizes  containing  (129  x  104),  (200  x  200)  and  (256  x  300)  points  in  the  radial  and 
azimuthal  directions  respectively,  to  establish  the  grid-independence  of  the  solution. 

Figures  (88)  through  (90)  display  several  flow  quantities  computed  on  the  different 
domains.  The  pressure  coefficient  and  vorticity  distribution  at  the  wall  (figures  88  and 
89),  which  together  determine  the  total  drag  on  the  cylinder,  differ  on  the  upstream 
side  of  the  cylinder  by  a  pointwise  maximum  of  2%  on  the  120 Rc  and  300 /?c  domains. 
Downstream  of  the  separation  point  ( 6  >  0.75 x),  on  the  portion  of  the  cylinder 
wall  facing  the  wake,  both  wall  pressure  coefficient  and  vorticity  are  independent  of 
computational  domain  size. 

The  streamwise  velocity  on  the  rear  flow  symmetry  axis  (figure  90),  is  influenced 
by  the  domain  size  in  a  region  extending  from  approximatly  10  diameters  downstream 
to  the  outer  domain  boundary.  In  the  near  wake  however,  the  velocity  distributions 
obtained  on  the  three  domains  are  identical.  The  pressure  coefficient  on  the  rear  flow 
axis,  shown  in  the  same  figure,  is  the  quantity  most  sensitive  to  domain  size.  On 
the  smaller  domain,  the  pressure  overshoots  and  increases  above  its  stagnation  value. 
As  the  domain  size  increases,  the  overshoot  diminishes  and  has  disappeared  on  the 
larger  domain.  Based  on  these  results,  the  domain  chosen  for  all  steady  and  unsteady 
laminar  computations  is  a  circle  of  radius  Rd  =  300i?c. 
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D.3  Laminar  Vortex  Shedding  at  Re  =  80, 100 

Two-dimensional  computations  at  Reynolds  numbers  80  and  100  are  performed  to 
establish  the  temporal  accuracy  of  the  numerical  method.  These  simulations  further 
test  the  far-held  boundary  condition  performance  in  unsteady  flow  by  examining 
whether  the  only  energetic  frequency  is  the  primary  wake  frequency.  Inaccuracies  in 
outflow  boundary  conditions  can  generate  modulating  frequencies  in  the  flow  response 
(Don  1989). 

Appendix  C  describes  the  family  of  grids  on  which  all  cylinder  computations 
are  performed.  Simulation  results  at  Reynolds  numbers  80  and  100  are  presented 
in  Sections  D.3.1  and  D.3. 2  respectively.  It  is  shown  that  the  Strouhal  frequency 
converges,  with  grid  refinement,  to  the  experimental  results  of  Williamson  (1991)  for 
two-dimensional  vortex  shedding.  The  performance  of  the  iterative  solver  is  evaluated 
in  Section  D.3.3  using  the  case  of  laminar  vortex  shedding  at  Reynolds  number  100. 
The  iterative  scheme  is  shown  to  have  good  convergence  properties  in  this  case. 

D.3.1  Vortex  Shedding  at  Reynolds  Number  80 

Table  (9)  in  Appendix  C  describes  the  grids  used  in  the  simulations  at  Reynolds 
number  80.  Radial  and  azimuthal  directions  are  refined  separately  in  this  study, 
while  the  computational  domain  radius  is  set  to  300 /?c.  Radial  grid  refinements  are 
performed  on  5  successive  meshes.  On  these,  the  radial  wall  spacing  decreases  from 
1.59 Rc  to  7.2  x  10~4RC,  the  stretching  factors  s\  and  s2  are  equal  and  fixed  at  1.03, 
and  the  65  total  azimuthal  points  are  equispaced.  The  angles  defining  both  end  points 
of  the  wake  envelop  are  equal  and  set  at  90  degrees,  resulting  in  standard  polar  mesh 
configurations. 

Azimuthal  refinements  of  the  wake  envelop  region  are  done  on  3  additional  grids, 
which  contain  respectively  48,  64  and  94  points  across  the  wake  layer.  The  last  grid 
features  a  wake  envelop  similar  to  that  presented  in  figure  (78).  It  further  uses  a 
variable  radial  grid  stretching  ranging  between  1.03  at  the  cylinder  surface  and  1.1 
at  the  outer  domain  boundary. 
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Grid  Number 

CFL  (Uoo) 

CFL  (Coo) 

UooAt/Rc 

St 

1 

1 

13 

0.254 

0. 

2 

6 

0.118 

0.137 

3 

12 

0.060 

0.148 

4 

39 

0.037 

0.149 

5 

225 

0.033 

0.150 

6 

178 

0.026 

0.152 

7 

154 

0.022 

0.152 

0.5 

72 

0.010 

8 

1 

154 

0.022 

Table  17:  Computed  Strouhal  number  at  Re=  80 

CFL  numbers  and  time-steps  used  on  each  grid  are  listed  in  table  (17).  The  time 
step  in  each  simulation  was  chosen  to  correspond  to  a  velocity-based  CFL  number 
of  1.  A  time-step  refinement  performed  on  grid  7  compares  solutions  computed  with 
CFL  numbers  of  1  and  0.5.  It  is  worthwhile  to  note  that  in  these  two-dimensional  sim¬ 
ulations,  the  implicit  time-marching  scheme  is  stable  and  accurate  for  CFL  numbers 
as  high  as  225. 

On  the  coarsest  mesh  (65  x  65),  no  vortex  shedding  can  be  sustained  (St  =  0).  As 
the  grid-density  increases  at  the  cylinder  surface,  the  Strouhal  number  reaches  a  value 
of  0.152.  Figure  94  displays  the  evolution  of  the  Strouhal  number  with  grid  refinement. 
It  indicates  that  the  value  St  =  0.152  is  closely  grid-independent.  This  result  is 
in  agreement  with  the  experimental  frequency  of  0.153  documented  by  Williamson 
(1991). 

On  the  densest  mesh  with  (320  x  127)  points,  a  reduction  of  the  time-step  by  a 
factor  of  2  does  not  affect  the  value  of  the  computed  Strouhal  number.  The  simulation 
on  the  mesh  (8)  with  the  wake  envelop  shown  in  figure  (78),  accurately  predicts 
the  Strouhal  frequency,  demonstrating  that  selectively  refining  the  wake  layer,  while 
maintaining  the  resolution  in  the  potential  region,  does  not  affect  the  accuracy  of  the 
vortex  shedding  frequency  at  this  Reynolds  number. 

The  least-squares  curve-fit  of  the  lift  coefficient  (Cl)  with  a  sine  function  is  shown 
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Grid  Number 

CFL  (U0 G) 

CFL  (Coo) 

Un&t/Rc 

St 

1 

1 

5 

0.0782 

0.153 

2 

6 

0.0458 

0.156 

3 

6 

0.0364 

0.157 

4 

5 

0.0262 

0.163 

5 

21 

0.0201 

0.164 

6 

61 

0.0159 

0.164 

Table  18:  Computed  Strouhal  number  at  Re=  100 

in  figure  (93).  The  maximum  pointwise  error  between  Cl  and  the  fitted  curve  is 
on  the  order  of  10-4  smaller  than  the  lift  amplitude,  indicating  that  no  energetic 
mode  other  than  the  primary  vortex  shedding  one  is  present.  Modulation  of  the  lift 
coefficient  response  by  frequencies  lower  than  the  Strouhal  frequency  is  a  numerical 
problem  which  can  arise  because  of  inaccuracies  in  far-field  boundary  conditions  (Don 
1989). 

Lift  and  drag  coefficient  decompositions  into  viscous  and  pressure  components  are 
shown  in  figures  (91)  and  (92).  At  this  Reynolds  number,  the  skin-friction  contributes 
appreciably  to  the  total  lift,  and  accounts  for  about  one-third  of  the  total  drag. 

D.3.2  Vortex  Shedding  at  Reynolds  Number  100 

The  grids  used  at  Reynolds  number  100  are  described  in  table  (10)  (Appendix  C). 
All  feature  an  azimuthal  wake  envelop  as  well  as  variable  radial  grid  stretching.  The 
cylinder  surface  and  outer  boundary  stretching  factors  are  fixed  at  Si  =  1.03  and 
s2  =  1.1  respectively.  Azimuthal  refinements  of  the  wake  envelop  are  performed 
on  grids  containing  40,  80,  100  and  130  points  across  the  wake  layer.  The  radial 
spacing  at  the  cylinder  surface  varies  on  the  6  grids  from  0.2078/4  to  0.0013/4.  The 
computational  domain  radius  is  fixed  in  all  cases  at  300 Rc. 

The  CFL  numbers  on  each  grid  are  summarized  in  table  (18).  Velocity-based  CFL 


Re 

st 

max  Cl„ 

max  Clv 

max  Cd„ 

max  Cdv 

80 

0.152 

0.215 

0.037 

0.976 

0.378 

100 

0.164 

0.297 

0.045 

1.021 

0.344 

Table  19:  Lift  and  drag  maxima 

numbers,  fixed  at  1  in  each  simulation,  define  the  associated  time-steps  and  sound- 
speed-based  CFL  numbers.  The  evolution  of  the  Strouhal  number  with  grid  refine¬ 
ment  is  displayed  in  figure  (98).  In  the  most  accurate  simulation,  the  Strouhal  fre¬ 
quency  stands  at  0.164,  in  agreement  with  the  experimental  value  of  0.164  (Williamson 
1991).  Lift  and  drag  coefficients  are  shown  in  figures  (95)  through  (97).  The  contri¬ 
bution  of  the  shear  stress  to  both  coefficients  (table  19)  is  smaller  at  Reynolds  number 
100  than  at  Reynolds  number  80.  The  maximum  lift  due  to  viscosity  is  13%  of  the 
total  lift  at  Reynolds  number  100,  and  15%  at  Reynolds  number  80.  The  maximum 
skin-friction  accounts  in  these  same  cases  for  respectively  25%  and  28%  of  the  total 
drag. 

Figure  (97)  demonstrates,  as  in  the  case  of  Reynolds  number  80,  the  fit  of  the  lift 
coefficient  response  with  a  sinusoidal  curve,  and  the  absence  of  energetic  low  frequency 
modes. 


D.3.3  Iterative  Solver  Performance 

The  performance  of  the  iterative  solver  is  tested  in  the  computation  of  the  Strouhal 
number  at  Reynolds  number  100.  The  fourth  grid  described  in  table  (10)  is  used  to 
demonstrate  the  behavior  of  the  residuals  with  two  different  time-steps.  In  both 
computations,  the  number  of  sub-iterations  per  time-step  is  set  to  3.  The  two  time- 
steps  correspond  to  velocity-based  Courant  numbers  of  1  and  0.5,  or  sound-speed 
Courant  numbers  of  5  and  2.5  respectively.  A  measure  of  convergence  error  in  the 
solution  at  time- level  n  and  sub-iteration  p  is  the  L\  norm  of  the  residual  error  || 
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(186) 


(equation  73,  Chapter  2),  defined  by: 

rLz  r27T  rRn 

IliT-i !  =  jf  J  J  |ir-p|ara^ 

The  decrease  in  this  residual  error  at  a  given  time  level,  expressed  by  the  ratio  of 
residual  norms  at  the  first  and  last  iterations  ||i?"-1||i/||i2n’3||i,  is  an  indicator  of  the 
performance  of  the  iterative  solver.  Time  evolutions  of  this  ratio  are  displayed  in 
figures  (99)  and  (100)  for  each  CFL  number.  Convergence  at  each  sub-iteration  is 
better  for  smaller  time-steps,  the  norm  of  the  residuals  dropping  faster  than  linearly 
with  decreasing  time-step.  At  CFL  number  of  5,  all  residuals  drop  by  a  factor  of 
about  1,000  in  three  iterations.  At  half  the  CFL  number,  the  residuals  drop  by  an 
additional  factor  of  approximately  4. 

The  iterative  method  with  three  sub-iterations  per  time-step  is  approximately 
three  times  faster  than  a  direct  inversion  method.  Because  of  the  modifications  made 
to  simplify  the  linearizing  Jacobians  (Chapter  2),  the  standard  proof  of  quadratic 
convergence  for  a  Newton  scheme  does  not  hold.  Empirically,  this  scheme  is  found  to 
be  robust.  The  accuracy  of  mean  surface  quantities,  including  Strouhal  number,  drag 
and  lift  coefficients,  does  not  seem  unduly  affected  by  fixing  the  number  of  iterations 
at  3,  incurring  at  each  time-step  a  finite  residual  error. 
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Figure  89:  Re=20;  Wall  vorticity 
-  :  Rd/Rc  =  60; - :  Rp/Rc  =  120; - :  Rd/Rc  =  300 
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Figure  91:  Re=80;  Lift  coefficient 
- :  Pressure  lift;  •  •  •  :  Viscous  lift;  -  :  Total  lift 
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Figure  92:  Re=80;  Drag  coefficient 
-  :  Form  drag;  •  •  •  :  Skin  friction;  -  :  Total  drag 
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Figure  93:  Re=80;  Strouhal  frequency  evaluation 
-  :  Computation;  •  :  Best  sinusoidal  fit 


Mesh  Size  xlu 


Figure  94:  Re=80;  Strouhal  frequency  convergence 
•  :  129  x  65;  A  :  200  x  65;  o  :  256  x  65 
x  :  320  x  65;  o  :  320  x  97;  +  :  320  x  127 


Appendix  E 

Linear  Stability  of  a  Forced 
Channel  Flow 


E.l  Formulation  of  the  Problem 


The  accuracy  of  the  numerical  treatment  of  the  convective  acceleration  is  tested  in 
a  periodic  channel  using  linear  stability  analysis.  The  work  necessary  to  drive  the 
fluid  through  the  channel  is  provided  by  an  external  periodic  force.  The  streamwise 
direction  in  the  channel  is  denoted  by  y,  and  x  is  the  cross-channel  direction. 

In  steady-state  incompressible  Poiseuille  flow,  the  maximum  streamwise  velocity 
is  denoted  v  =  {/<»,  and  is  related  to  the  the  pressure  gradient  across  the  channel  by: 

I  -  -h*”  (187) 


where  h  is  the  channel  half-width  and  p,  the  viscosity.  The  quantity 


1  dp 
Poo  dy 


2 

pooh 2 


liUoo 


(188) 


is  a  constant  streamwise  force  per  unit  mass,  and  represents  the  source  of  momentum 
necessary  for  the  flow  to  achieve  the  steady-state  velocity  r  =  !/«,  on  the  channel 
centerline. 

I 
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The  thermodynamic  pressure  distribution  obtained  by  integrating  (188)  is  linear 
in  y.  It  is  thus  not  a  solution  to  the  boundary- value  problem  describing  the  com¬ 
pressible  flow  in  a  periodic  channel.  To  circumvent  this  difficulty,  an  external  forcing 
is  introduced  to  drive  the  flow,  in  the  form  of  a  momentum  and  energy  flux  calibrated 
to  match  the  magnitude  of  the  force  given  in  (188).  The  external  force  imposed  in 
the  streamwise  momentum  equation  is  a  constant  denoted  by  A.  Corresponding  to 
this  momentum  source,  an  energy  source,  equal  in  magnitude  to  the  product  of  A 
and  the  streamwise  velocity  v,  appears  in  the  energy  equation.  The  flux  that  must 
be  imposed  as  an  external  source  term  in  the  Navier-Stokes  equations  written  in 
conservative  form  is  then 


T  =  A  1 


(189) 


The  constant  A  is  defined  by: 


2JMc 

Re 


(190) 


where  the  reference  velocity  is  the  sound  speed  c^.  The  symbol  J  refers  to  the 
Jacobian  of  the  transformation  from  physical  to  computational  space  (equation  42, 
Chapter  2).  Upon  introduction  of  the  external  forcing,  the  governing  equations  (equa¬ 
tion  45,  Chapter  2)  become: 


—6  +  —Ft  +  —F  +  —F  =  f 

dtQ+dC*+dri  v  +  dz  z 


(191) 


The  vector  Q  represents  the  conservative  variables  (equation  44,  Chapter  2).  Its 
definition  is  restated  here  for  clarity: 


Q  =  93  =  J\  pv 


(192) 
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It  is  worth  noting  that  although  the  external  force  A  itself  is  known  and  constant, 
the  external  energy  source  involves  the  streamwise  component  of  momentum  ( v ). 
The  time  integration  of  the  equations  above  thus  involves  both  explicit  and  implicit 
contributions  by  the  external  forcing.  The  corresponding  implicit  Jacobian  of  the 
forcing  term  is  described  in  the  next  Section. 


E.2  Numerical  Implementation 

The  modifications  to  the  time-marching  schemes  arising  because  of  the  forcing  T  are 
derived  by  considering  the  time-marching  scheme  given  by  (73)  in  Chapter  2: 

AQn+1  -  c3^A Qn+1  =  (c,  -  1  )Qn  +  c2<5"-1  +  c3jtQn  (193) 

where  n  represents  the  time  level,  A  Qn+1  =  Qn+1-Qn,  and  the  coefficients  ci  through 
c3  are  functions  of  the  time  step.  Adding  the  forcing  flux  T  to  the  right-hand  side  of 
the  Navier-Stokes  equations  results  in  equation  (191),  which  substituted  into  (193) 
leads  to: 

A<T+1  +  C3A(^  +  +  Th*’)  " 

implicit  contribution 

(c1-1)0"  +  cj(3’'-1-C3(^F(  +  ^F„  +  £f.)"+  gP  (194) 

explicit  contribution 

The  implicit  contribution  Af  of  the  forcing  flux  is  a  new  term  which  must  be  expressed 
as  a  function  of  the  solution  vector  A Q.  It  is  linearized  in  time  using  a  Taylor  series 
expansion: 

Afn+1  =  ^)nAgn+1  +  i?(Afn+1)2  (195) 

where  the  flux  Jacobian  is  obtained  directly  from  the  definition  of  T(Q)  given  in  (189) 


t 
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0  0  0  0  0 

0  0  0  0  0 

=  A  0  0  0  0  0 

0  0  0  0  0 


(196) 


\-q3/qj  0  1/ft  0  0/ 

This  expression  for  the  Jacobian  of  the  forcing  completes  the  definition  of  the 
additional  terms  appearing  in  the  time-marching  scheme.  The  non-linear  system  of 
equations  is  inverted  using  the  iterative  technique  outlined  in  Chapter  2. 

E.3  Steady  One-Dimensional  Solution 


The  governing  steady-state  one-dimensional  compressible  equations  in  physical  space 
are  obtained  from  (191)  as: 

+  =  f  (197) 


with  J  =  1  and  u(ar)  =  w(x )  =  0.  In  terms  of  the  primitive  variables,  the  equations 


reduce  to: 


(7  -  l)(e  -  pv2/ 2) 

—vx/Re 

ire{vvx  +  lTx/Pr(i  -  1)) 


=  A  1 


(198) 


Channel  walls  are  located  at  x  =  ±1,  where  the  boundary  conditions  are: 


u  =  u  =  uj  =  0,  T  —  Tw  =  — . 

7 


(199) 


Substituting  the  definition  of  A  (equation  190),  the  solution  of  this  boundary-value 
problem  is: 

v(x)  —  Moo(l  —  x 2)  (200) 

,  .  3  /AA.  \ 


'  PrMK 7  -  1)(1  -  x4)  +  3 7rt 


(201) 


(202) 
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In  this  solution,  the  pressure  is  constant,  p  =  I/7,  the  streamwise  velocity  is  formally 
identical  to  the  incompressible  Poiseuille  solution,  and  the  temperature  is 

T(x)  =  =  Tw  +  —Ml(i  -  1)(1  -  x4).  (203) 

7 P\x)  J7 

In  the  limiting  case  of  zero  Mach  number,  the  density  becomes  constant  and  the 
solution  is  well  behaved.  The  governing  equations  generate  this  solution  because  the 
molecular  viscosity  is  assumed  independent  of  temperature,  and  because  the  forcing 
(A)  is  per  unit  volume  instead  of  per  unit  mass. 


E.4  Small-Disturbance  Equations 


Substituting  perturbed  primitive  variables  into  the  conservative  form  of  the  Navier- 
Stokes  equations  (191)  written  in  cartesian  coordinates,  and  neglecting  all  products 
of  perturbation  fields,  leads  to  the  three-dimensional  small-disturbance  equations, 


d_ 

dt 


P'  \ 
pu' 

pv'  +  p'v 
pw' 
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V  e'  / 


f  pu'  \ 
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J_d_ 
Re  dz 


K  +  < 

w'v  +  v'x 


(7  —  1  )RePr 


=  A  0  I  . 


-§«  +  <) +  K 

\  v(w'y  +  v'z )  /  \1 

where  the  perturbation  pressure  and  temperature  are 


T'  +T'  +  T 

x  xx  1  yy  '  zz 


1  /  /  1  /  2 

•p  —  e  —  pv  v  —  -p  v 


(204) 


P\  ~iPJ 


At  the  solid  wall,  we  impose 


v!  —  v'  =  w'  =  0,  T'  =  0. 


(205) 


(206) 


These  are  the  only  physical  boundary  conditions  required,  since  the  normal  momen¬ 
tum  equation  at  the  wall  defines  the  normal  gradient  of  pressure  at  the  wall.  The 
channel  geometry  is  periodic  in  the  streamwise  ( y )  and  spanwise  ( z )  directions,  and 
disturbances  are  decomposed  in  Fourier  modes  along  these  directions.  Considering 
a  generic  Fourier  mode,  we  seek  solutions  to  the  small-disturbance  equations  of  the 


v'  ( x,y,z,t ) 


jaiy+pz-ct) 


(207) 


Substitution  into  the  small-disturbance  equations  yields  an  eigenvalue-eigenvector 
problem  for  c, 

L2 ^  +  L +  Lo$  =  c$  (208) 

dx 2  Ox 
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where  \P  is  the  solution  eigenvector 


/u\ 

A 

V 

A 

w 

P 

\T ) 


(209) 


The  matrix  L0  is  proportional  to  the  product  of  two  matrices  Mi  and  M2 


Lq  — - M\  M2 

ap 


where  Mi  and  M2  are  given  by 
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(210) 
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Matrices  L\  and  L2  are  given  by 

(211) 


(212) 


For  given  wave-numbers  a  and  /?,  the  solution  of  this  eigenvalue  problem  yields  the 
eigenfunctions  (u,v,w,p,T)  and  the  complex  frequency  c.  The  perturbation  energy 
is  defined  as 

E(t)  =  £Z  £*  f  ( u ,2  +  i>'2  +  w'2)dxdydz.  (213) 

Its  evolution  in  time  is  exponential, 

E(t)  =  E(0)e2cit  (214) 

where  Cj  is  the  imaginary  part  of  the  eigenvalue  c.  The  objective  of  the  linear-stability 
simulations  is  to  test  the  code  accuracy  in  predicting  the  energy  growth-rate  (2c,). 


/  0 


Xi  = 


1 

3  pRe 

_g_ 

3  pRe 


ihzH 

<*P1 


1 

3  pRe 

o 

o 


..  A .  _L  n\ 

3  pRe  ap  1 


-M(T  -  life 
-2<(7  -  D* 


and 


Lo  = 


aRe 


0  0 
0  0 
0  0 
0  0/ 


/4/3  p  0  0  0  0  \ 

0  1/p  0  0  0 

0  0  1/p  0  0 

0  0  0  0  7 /Pr 

V  0  0  0  0  7  j  pPr  / 


E.5  Three-Dimensional  Eigenfunctions 

The  oblique  mode  a  =  /3  =  1  is  selected  for  three-dimensional  simulations.  The 
equations  do  not  admit  growing  modes  in  this  case,  and  the  slowest  decaying  Orr- 
Sommerfeld  mode  is  chosen  as  the  perturbation.  The  corresponding  frequency  is 
c  =  0.02962395  —  iO. 00220154  at  Reynolds  number  7,500  and  Mach  number  0.1. 
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Scaling  Factor 

Real  Part 

Imaginary  Part 

m 

0.08143 

0.07145 

0.9964 

0.5974 

l^lmax 

0.81796 

0.7994 

|F|max 

0.00274 

0.00283 

|T|max 

0.01504 

0.01606 

Table  20:  Scaling  factors  of  3-D  eigenvalues 


The  corresponding  incompressible  solution  gives  as  the  slowest  decaying  mode  c  — 
0.02963201  —  *0.00219390,  a  difference  in  norm  of  0.025%. 

The  real  and  imaginary  parts  of  the  eigenfunctions  are  shown  in  figures  (101)  and 
(102).  These  functions  are  scaled  on  their  respective  maxima  across  the  channel, 
which  are  listed  in  table  (20).  The  results  of  the  simulations  performed  are  described 
in  Chapter  2.  In  all  simulations,  the  point  distribution  across  the  channel  is  given  by 
a  hyperbolic  tangent  with  a  stretching  factor  a: 


tanh(o:x?) 

tanh(o) 
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j  Re  $ 


Figure  103:  Numerical  dissipation  in  steady  channel  flow  solution 
Stretching  factors: -  :  a  =  2; - :  a  =  3;  •  •  •  :  a  =  4.5 
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Appendix  F 

Numerical  Aspects  of  LES 


F.l  Introduction 

Large-eddy  simulations  are  based  on  the  division  of  turbulent  flows  into  large  and 
small-scale  motions.  The  energetics  of  the  large  scales  are  computed  directly,  whereas 
the  effect  of  the  small  scales  on  the  large  scales  is  modeled.  The  performance  of  this 
division  of  labor  between  the  filtered,  large-scale  equations  and  the  subgrid-scale 
model  hinges  partially  on  the  ability  of  the  mesh,  given  a  numerical  method,  to 
accurately  resolve  the  dynamically  important  motion.  This  appendix  summarizes  the 
large-eddy  grid-refinement  study  performed  at  a  Reynolds  number  of  3,900  based  on 
free-stream  velocity  and  cylinder  diameter.  The  objective  is  to  select  a  grid  on  which 
comparisons  between  simulations  using  different  subgrid-scale  eddy- viscosity  models 
can  be  established.  All  calculations  presented  in  this  appendix  were  performed  with 
the  fifth-order  accurate,  upwind-biased  scheme  for  the  convective  derivatives. 

The  subgrid-scale  eddy  viscosity,  given  by  (30)  in  Chapter  2,  is  of  the  form  ut  — 
CA2p  |5*|,  where  C  is  a  function  of  space  depending  on  the  Leonard  stress  tensor, 
density  and  rate  of  strain  fields.  The  convergence  properties  of  the  subgrid-scale 
Reynolds  stresses  as  the  grid  is  refined  have  not  been  thoroughly  investigated.  The 
purpose  of  the  present  study  is  thus  to  examine  the  variations  of  the  mean  quantities 
in  the  wake  and  at  the  cylinder  surface  as  grid  resolution  increases.  Most  velocity 
and  resolvable  Reynolds  shear  stress  statistics  near  the  wake  centerline  are  found 
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to  oscillate  with  grid  refinement,  indicating  that  the  finer  meshes  capture  the  main 
statistical  flow  features. 

The  following  section  discusses  three-dimensional  initial  conditions.  Grid  and 
time-step  refinement  studies  are  presented  in  section  F.3.  The  final  mesh  selected  for 
large-eddy  simulations  is  described  in  section  F.4. 


F.2  Three-Dimensional  Initial  Conditions 

Three-dimensional  calculations  are  initialized  with  planes  of  data  obtained  from  pre¬ 
liminary  two-dimensional  simulations,  each  plane  corresponding  to  the  flow-field  at  a 
different  phase  of  the  vortex  shedding  cycle.  The  resulting  initial  fields  numerically 
satisfy  the  continuity  equation,  and  the  three-dimensional  motions  are  generated  by 
a  spanwise  force  imbalance.  This  technique  is  used  to  generate  turbulent  fields  only 
once  for  each  spanwise  box-size.  Simulations  performed  on  computational  domains 
having  identical  spanwise  lengths  but  different  grids  are  initialized  with  fields  inter¬ 
polated  from  solutions  computed  on  coarser  meshes  to  minimize  the  perturbation  of 
the  large  scales. 

In  the  present  simulations,  no  external  perturbations  are  imposed  upstream  of  the 
cylinder.  In  the  absence  of  external  perturbations,  the  three-dimensionality  of  the  flow 
was  found  to  be  numerically  self-sustained  for  spanwise  grid  spacings  smaller  than 
Az  =  QARC,  corresponding  to  5  spanwise  points  per  cylinder  diameter.  For  coarser 
spanwise  resolutions  the  flow  reverts  to  a  two-dimensional  state  after  a  transient  of  1 
to  2  vortex  shedding  cycles. 
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Grid  size  (r,  6 ,  z) 

U^At/Rc 

CFL(Uoo) 

CFL(coo) 

88  x  90  x  32 

0.03 

1.2 

15 

88  x  90  x  48 

0.01 

0.4 

5 

116  x  136  x  48 

0.005 

0.3 

6 

174  x  128  x  48 

0.005 

0.7 

13 

Table  21:  Time  integration  parameters 

F.3  Spatial  and  Time  Resolution  Study 

F.3.1  Grid  Sizes 

The  selection  of  a  box-size  in  the  spanwise  direction  (z)  is  such  that  the  two-point 
correlation  function  of  velocity  component  u, 

l  ^ 

Ru(0,0,r)  =  m(z,y,z,t)  m(x,y,z  +  r,t)  -  u\D{x,y,t)dt  -  u]{x,y)  (217) 

is  close  to  zero  at  the  box  center.  The  bar  over  a  quantity  in  that  expression  de¬ 
notes  an  average  over  time  and  spanwise  distance.  The  two-dimensional  component 
of  the  velocity  Ui2D  is  by  definition  independent  of  the  spanwise  direction  (chapter  2, 
section  2.8).  This  relation  assumes  that  the  random  and  two-dimensional  compo¬ 
nents  of  the  flow  are  uncorrelated  (Cantwell  k  Coles  1983).  Recent  experiments 
have  investigated  the  size  and  nature  of  the  structures  present  in  the  intermediate- 
wake  behind  the  cylinder.  These  large-scale  three-dimensional  structures  are  pairs  of 
counter-rotating  streamwise  vortices  (Hayakawa  k  Hussain  1989).  Between  Reynolds 
numbers  330  and  21,000,  the  mean  spanwise  spacing  of  these  structures  is  constant 
at  about  one  pair  per  diameter  (Bays-Muchmore  k  Ahmed  1993).  The  correlation 
functions  R„(0, 0,r)  were  however  not  documented  in  these  experiments. 

The  only  location  where  experimental  measurements  are  available  to  evaluate 
spanwise  two-point  correlations  is  three  diameters  downstream  of  the  cylinder.  Ong 
k  Wallace  (1994)  computed  the  spanwise  two-point  correlation  of  the  hot-wire  voltage 
E  at  that  location  at  Reynolds  number  3, 900,  using  a  single-wire  probe.  The  effects 
of  both  streamwise  and  vertical  velocities  are  thus  included  in  these  correlations.  The 
measurements,  taken  on  the  wake  centerline  (y  —  0),  are  shown  in  figure  104.  They 

j 

i 
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represent  the  sum  of  the  correlation  function  and  the  squared  periodic  component  of 
the  voltage  fluctuation: 

R*ee(  0, 0,  r)  =  Ree(  0, 0,  r)  +  i  £  E22D(x,  y,  t)dt  (218) 

The  two-dimensional  component  of  the  voltage  fluctuations  is  not  available  experi¬ 
mentally.  To  derive  the  most  pessimistic  estimate  of  the  location  at  which  the  correla¬ 
tion  Ree  is  zero,  we  assume  that  the  voltage  function  E(t )  is  defined  uniquely  by  the 
streamwise  velocity  u(t),  and  that  the  periodic  component  is  the  net  two-dimensional 
component  of  the  flow.  At  x/D  =  3,  the  experimental  value  of  U2l>/#*u(0,0,0)  is 
about  0.1.  The  correlation  function  i?uu(0,0,r)  at  three  diameters  downstream  is 
thus  zero  at  approximately  r/D  =  1.5.  In  all  the  simulations  documented  in  this 
report,  the  spanwise  box-size  is  set  to  twice  this  distance,  at  Lz  =  irD. 

The  layouts  of  the  meshes  used  in  this  resolution  study  are  discussed  in  section  C.3 
of  Appendix  C.  To  establish  that  the  final  mesh  selected  resolves  the  main  features  of 
the  flow,  mean  velocities  and  Reynolds  shear  stresses  are  presented  below  on  4  grids 
containing  (88  x  90  x  32),  (88  x  90  x  48),  (116  x  136  x  48)  and  (174  x  128  x  48)  points 
in  the  radial,  azimuthal  and  spanwise  directions  respectively. 

F.3.2  Time  Stepping  and  Iteration  Residuals 

The  time-steps  and  CFL  numbers  corresponding  to  each  calculation  are  displayed  in 
table  (21).  Statistics  are  compiled  in  each  simulation  over  60  to  65  time  units  based 
on  cylinder  radius  and  free-stream  velocity,  corresponding  to  approximately  6  vortex 
shedding  cycles.  The  convergence  of  the  iteration  residuals,  illustrated  in  the  case 
of  the  (88  x  90  x  48)  simulation,  is  shown  in  figure  (105).  This  figure  displays  time 
histories  of  the  residual  errors  in  each  governing  equation,  defined  in  section  D.3  of 
Appendix  D,  over  approximately  half  a  vortex  shedding  cycle.  The  number  of  sub¬ 
iterations  per  time-step  is  fixed  at  3.  Residuals  drop  on  the  average  by  a  factor  of  103 
at  each  time-step,  with  extrema  in  all  five  equations  of  motion  as  low  20  and  as  high 
as  104.  Although  convergence  histories  at  Reynolds  number  3,900  are  noisier  than 
those  at  Reynolds  number  100  (figures  99  and  100,  Appendix  D),  the  level  of  mean 
convergence  is  comparable  in  both  cases. 
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F.3.3  Mean  Convergence  with  Time  of  Integration 

The  convergence  of  mean  flow  statistics  with  time  of  integration  is  examined  in  the 
case  of  the  (116  x  136  x  48)  simulation.  Table  (23)  summarizes  the  evolution  of  the 
main  near-wake  features  integrated  over  6  vortex  shedding  cycles.  Figures  (106)  and 

(107)  display  the  wall  pressure  coefficient  and  shear  stress,  as  well  as  velocity  profiles 
inside  and  downstream  of  the  recirculation  bubble  at  different  times  of  integration. 
Total  drag,  skin-friction,  upstream  separation  angle  and  back-pressure  are  close  to 
their  final  values  after  two  shedding  cycles.  Other  quantities,  including  secondary 
separation  angles,  bubble  length,  wall  stress  and  wake  velocities  require  up  to  5  shed¬ 
ding  cycles  for  convergence,  corresponding  to  integrations  over  approximately  50i?c/ 
Uoo  time  units. 

F.3.4  Effect  of  Refinement  on  Mean  Wake  Quantities 

Statistics  obtained  from  each  simulation  are  summarized  in  table  (24)  and  figures 

(108)  through  (116).  Table  (24)  lists  mean  quantities  of  interest  at  the  cylinder 
surface  and  in  the  main  recirculation  bubble.  These  parameters  are  observed  to 
either  oscillate,  or  decrease,  reach  a  minimum  and  increase,  or  undergo  the  reverse 
trend,  as  the  grid  is  refined.  In  particular,  the  oscillatory  behavior  of  the  back¬ 
pressure  coefficient  and  separation  angles  can  be  observed  in  figure  (108),  which 
displays  pressure  coefficient  and  vorticity  around  the  cylinder  circumference  for  each 
simulation. 

Three  of  the  most  important  wall  quantities  are  the  back-pressure  coefficient,  the 
upstream  separation  location  of  the  boundary  layer  and  the  total  drag  coefficient. 
These  all  lay  within  experimental  error  in  the  two  most  resolved  simulations.  The 
back-pressure  coefficients  are  equal  and  stand  at  Cpb  =  —0.94  in  these  simulations, 
while  the  upstream  separation  angle,  equal  to  6\  =  86.2°  in  the  most  accurate  case, 
differs  by  about  1%  from  is  value  on  the  116  x  136  x  48  mesh.  The  experimental 
value  of  the  drag  coefficient  is  Cp  =  0.98  i  0.05  at  Reynolds  number  3,900.  The  most 
accurate  simulation  predicts  a  drag  coefficient  of  Cd  =  0.98.  The  Strouhal  frequency 
however  is  underestimated  on  all  meshes,  and  stands  on  the  finest  one  at  0.20,  5% 
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below  the  experimental  bound  of  0.215  ±  0.003  (Cardell  1993).  On  the  upstream  face 
of  the  cylinder,  the  laminar  boundary-layer  velocity  profiles  ve(r)  are  shown  in  figure 
(109)  at  angles  6  =  90,120  and  150  degrees.  Although  the  boundary-layer  velocities 
are  underestimated  on  the  mesh  with  32  spanwise  points,  they  are  well  resolved  in 
the  48-spanwise-point  simulations. 

One  of  the  most  significant  near-wake  statistics  is  the  length  of  the  recirculation 
bubble  downstream  of  the  cylinder.  It  equals  1.18  ±  0.05 D  according  to  Lourenco’s 
PIV  experiment  (1993),  and  1.33  ±  0.2  based  on  Cardell’s  hot-wire  investigation 
(1993).  Its  computed  values,  illustrated  by  bubble- velocity  and  pressure  distribu¬ 
tions  in  figure  (110),  in  the  two  most  accurate  simulations  are  1.36Z?  and  1.44D, 
which  is  within  Cardell’s  error  range  but  higher  than  the  PIV  estimates.  Mean  ve¬ 
locity  and  Reynolds  shear  stress  profiles  across  the  wake  are  displayed  in  figures 
(111)  through  (116)  at  several  locations  downstream.  Reynolds  shear  stress  and  ve¬ 
locity  local  extrema,  oscillate  as  resolution  increases.  Streamwise  velocities  exhibit 
±5%  oscillations  of  their  maximum  norms  based  on  local  centerline  velocity  deficit. 
Peak  values  of  Reynolds  shear  stress  vary  in  a  10%  band  in  the  (88  x  90  x  48)  and 
(174  x  128  x  48)-simulations. 

F.4  Final  Grid  Selected 

These  results  indicate  that  on  the  (116  x  136  x  48)-point  mesh,  near- wake  mean  ve¬ 
locities  and  Reynolds  stresses  are  resolved  within  uncertainty  bounds  of  5%  and  10% 
respectively.  This  mesh  is  thus  appropriate  for  turbulence  simulations.  However  in 
two  computations  with  and  without  the  dynamic  subgrid  scale  model,  differences  in 
mean  velocities  were  small  and  confined  to  the  first  2  diameters  of  the  wake.  To  en¬ 
hance  the  impact  of  the  subgrid  scale  model,  the  final  mesh  selected  is  similar  to  the 
(116  x  136  x  48)  grid  with  three  modifications:  (1)  the  radial  resolution  at  the  cylinder 
surface  is  chosen  to  match  that  of  the  (174  x  128  x  48)  mesh;  (2)  the  resolved  wake 
region  is  extended  from  5  to  10  diameters  downstream;  (3)  downstream  of  the  bubble 
closure  point,  the  mesh  is  stretched  geometrically  in  the  radial  direction.  The  next 
section  demonstrates  that  this  stretching  has  a  minimal  effect  on  the  mean  velocity, 
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while  the  eddy  viscosity  increases  locally  by  up  to  40  percent.  These  three  modifica¬ 
tions  translate  into  a  mesh  containing  (144  x  136  x  48)  points  in  the  radial,  azimuthal 
and  spanwise  directions  respectively.  Its  physical  characteristics  are  documented  in 
Appendix  C.  The  grid  layout  is  displayed  in  figure  (86). 

F.4.1  Impact  of  Radial  Grid  Stretching  Past  the  Bubble 
Closure 

In  the  resolved  wake  downstream  of  the  recirculation  bubble  (L  <  r  <  Rw ) ,  meshes 
used  in  the  refinement  study  above  have  a  constant  radial  spacing,  (figure  80).  In 
the  final  grid  selected  described  above,  the  mesh  stretches  radially  downstream  of  the 
recirculation  bubble  in  order  to  enhance  the  contribution  of  the  turbulence  model 
within  the  resolved  wake  region.  The  objective  of  this  section  is  to  demonstrate  that 
this  radial  coarsening  of  the  mesh  increases  the  turbulent  eddy  viscosity  relative  to 
the  molecular  viscosity  without  unduly  affecting  the  accuracy  of  the  solution. 

The  impact  of  radial  grid  stretching  for  r  <  Rw  is  demonstrated  using  the  (116  x 
136  x  48)  simulation  described  in  tables  (11)  and  (12)  as  a  base  case.  Geometric 
stretching  in  the  radial  direction  between  1  and  5  diameters  downstream  with  a 
factor  of  1.02  results  in  a  (106  x  136  x  48)-point  mesh.  Stretched  and  uniform  radial 
grid-spacing  distributions  are  displayed  in  figure  (85).  Azimuthal  and  spanwise  grid 
construction  parameters  are  identical  to  those  of  the  unstretched  (116  x  136  x  48)- 
point  mesh.  On  both  grids,  the  initial  condition  is  a  fully  developed  field.  Statistics 
are  averaged  over  l\ARc/Uoo  time  units,  corresponding  to  approximately  1  vortex 
shedding  cycle.  The  time-step  is  0.005i?c/f/oo,  which  translates  into  velocity  and 
sound-speed  based  CFL  numbers  of  0.3  and  6  respectively. 

Table  (22)  lists  mean  quantities  of  interest  at  the  cylinder  surface  computed  on 
both  grids.  Total  drag  (Cd  =  0.98)  and  skin-friction  (C /  =  0.90  x  10-2)  are  identical 
in  both  cases.  Differences  in  back-pressure  coefficient  and  upstream  separation  angles 
stand  at  1%  and  0.1%  respectively. 

Mean  eddy  viscosity  and  velocities  at  various  locations  in  the  wake  are  displayed 
in  figures  (117)  and  (118).  The  maximum  difference  in  mean  velocity  occurs  3.5 
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Wake  mesh 

cPb 

CD 

Cf  x  100 

0i 

Unstretched 

-0.90 

0.98 

0.90 

84.6 

Stretched 

-0.91 

0.98 

0.90 

84.5 

Table  22:  Aerodynamic  coefficients  on  stretched  and  unstreched  meshes 

diameters  downstream,  where  the  streamwise  velocity  difference  is  approximately  3% 
of  the  centerline  velocity  deficit.  On  the  stretched  mesh,  the  eddy  viscosity  peaks  up 
to  40%  higher  than  on  the  uniform  wake  grid,  indicating  that  the  radial  stretching 
and  coarsening  of  the  grid  have  the  targeted  property  of  enhancing  subgrid-scale 
dissipation  while  having  little  effect  on  the  mean  flow. 
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Table  23:  Sampling  error  analysis;  (116  x  136  x  48)  simulation 
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Figure  106:  Wall  pressure  and  vorticity  convergence 
Shedding  cycles:  •  :  1; - :  2;  •  •  •  :  3; - :  4; - :  5; 
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Table  24:  Near- wake  result  summary 
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Figure  108:  Wall  pressure  and  vorticity 

- :  88  x  90  x  32;  •  •  •  :  88  x  90  x  48 

-  :  116  x  136  x  48;  -  :  174  x  128  x  48 


Figure  109:  Tangential  velocity  radial  profiles  at  the  cylinder  surface 

- :  88  x  90  x  32;  •  •  •  :  88  x  90  x  48; - :  116  x  136  x  48 

-  :  174  x  128  x  48  (a):  6  =  90°;  (b) :  6=  120°;  (c) :  6  =  150° 
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Figure  111:  Streamwise  velocity  u/f/oo  at  x/D  =  0.58  (a),  1.06  (b),  1.54  (c),  2.02  (d) 

- .  88  x  90  x  32;  •  •  •  :  88  x  90  x  48; - :  116  x  136  x  48 

-  :  174  x  128  x  48 


Figure  112:  Streamwise  velocity  u/Uoo  at  x/D  =  2.50  (a),  3.06  (b),  3.54  (c),  4.50  (d) 

- :  88  x  90  x  32;  •  •  •  :  88  x  90  x  48; -  :  116  x  136  x  48 

-  :  174  x  128  x  48 


Figure  113:  Vertical  velocity  v/Uc <>  at  x/D  =  0.58  (a),  1.06  (b),  1.54  (c),  2.02  (d) 


Figure  114:  Vertical  velocity  v/Uoo  at  x/D  =  2.50  (a),  3.06  (b),  3.54  (c),  4.50  (d) 

- :  88  x  90  x  32;  •  •  •  :  88  x  90  x  48 

- :  116  x  136  x  48;  -  :  174  x  128  x  48 


Figure  115:  Total  Reynolds  shear  stress  in  the  formation  zone 
x/D  =  0.58  (a),  1.06  (b),  1.54  (c),  2.02  (d) 

- :  88  x  90  x  32;  •  •  •  :  88  x  90  x  48; - :  116  x  136  x  48 

-  :  174  x  128  x  48 
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